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Abstract 

Small,  lightweight  aircraft  are  of  increasing  interest  to  the  Air  Force,  for  addressing  poten¬ 
tial  threats  in  urban  environments  and  other  complex  terrain.  In  order  to  avoid  obstacles, 
respond  to  gusts,  and  track  potentially  elusive  targets,  highly  agile  maneuvers  will  be  re¬ 
quired,  for  which  the  standard  quasi-steady  aerodynamic  models  are  not  accurate.  This 
project  addresses  the  development  of  a  hierarchy  of  models  for  unsteady  aerodynamics,  in 
a  framework  that  is  suitable  for  control  design.  The  ultimate  goal  is  to  use  these  models  for 
the  design  of  flight  controllers,  for  instance  to  accurately  track  a  trajectory  in  the  presence 
of  large  disturbances. 

Classical  aerodynamic  models  by  Theodorsen  and  Wagner  were  shown  to  be  accurate  in 
many  cases,  but  break  down  at  high  angles  of  attack  or  for  very  rapid  maneuvers.  Finite- 
time  Lyapunov  exponents  (FTLE)  were  used  to  identify  meaningful  flow  structures  and 
elucidate  the  flow  physics.  Methods  were  developed  to  perform  balanced  model  reduction 
without  the  need  for  adjoint  information,  and  for  systems  with  periodic  motions  (such  as 
vortex  shedding).  These  were  used  to  develop  models  of  unsteady  pitching  and  plunging 
maneuvers  that  are  well-suited  to  control  design  techniques.  The  methods  were  developed 
for  linearized  models  (small-amplitude  motions)  but  also  perform  well  for  large-amplitude 
maneuvers. 

Contents 


1  Executive  summary 

1.1  Techniques! . 

1.2  Summary  of  this  report 

2 

Background:  classical  models 

2.1 

Theodorsen’s  frequency  response . 

2.2 

Wagner’s  indicial  response 

23 

23 

26 


2 


15  Balanced  model  reduction!  30 

5.1  Balanced  POD  models  without  adjoints| . 31 

5.2  Balanced  truncation  for  periodic  systems . 36 


6 

Control-oriented  extensions  of  classical  unsteady  models 

6.1  Reduced-order  approximations  of  Wagner’s  indicial  response 

6.2 

Results:  small- amplitude  maneuvers 

. 54 

6.3 

Results:  large-amplitude  maneuvers 

. 56 

16.4  Conclusions! . 59 


1  Executive  summary 

Most  models  used  for  flight  control  in  today’s  aircraft  assume  that  the  forces  and  moments 
on  the  aircraft  are  quasi-steady:  they  depend  only  on  the  velocity  of  the  vehicle  relative 
to  the  surrounding  air.  However,  unsteady  effects  become  increasingly  important  for  Micro 
Air  Vehicles,  where  rapidly  changing  gusts  and  large  accelerations  render  the  quasi-steady 
assumption  invalid. 

The  purpose  of  this  project  is  to  obtain  improved  models  for  unsteady  aerodynamics,  in 
a  form  that  is  suitable  for  control  design.  Classical  approaches  to  unsteady  aerodynamics 
include  the  frequency-domain  models  of  Theodorsen  [58]  and  the  time-domain  models  of  [63] , 
but  these  are  not  directly  suitable  for  control  design.  A  handful  of  more  modern  studies 
are  available  for  rolling  delta  wings  [4,  34]  and  for  airfoils  undergoing  dynamic  stall  [18L 141], 
but  these  still  have  many  limitations.  For  instance,  the  model  in  |18|  must  be  calibrated 
to  experimental  data,  and  often  does  not  match  data  it  was  not  specifically  calibrated 
against  [9].  The  model  m  is  linear,  and  any  nonlinear  effects  such  as  bifurcations  and 
hysteresis  are  not  captured.  The  focus  of  the  present  effort  is  to  obtain  models  that  overcome 
these  various  shortcomings,  and  may  be  used  for  designing  flight  controllers. 

1.1  Techniques 

The  overall  goal  of  this  work  is  to  provide  improved  models  for  design  of  flight  controllers. 
Our  focus,  however,  is  on  the  aerodynamic  models:  predicting  lift,  drag,  and  moments  from 
quantities  such  as  freestream  velocity  and  angle  of  attack.  In  traditional  flight  dynamic 
models,  these  lift  and  drag  forces  are  typically  treated  as  quasi-steady:  for  instance,  lift 
is  treated  as  a  static  function  of  angle  of  attack  (equal  to  the  steady-state  value).  When 
unsteady  aerodynamics  become  important,  this  approach  is  not  sufficient,  and  so  we  treat 
the  aerodynamic  models  separately,  as  shown  in  Figure  [l]  The  focus  of  the  present  effort  is 
to  determine  models  for  the  “Aerodynamics”  block  in  Figure  [T] 

The  models  we  are  developing  in  this  project  build  on  a  number  of  numerical  tools,  which 
we  describe  briefly  here. 

Immersed  boundary  solver.  The  starting  point  for  many  of  the  modeling  approaches  we 
employ  is  a  high-fidelity  direct  numerical  simulation  (DNS).  (In  many  cases,  experimental 
data  is  equally  suitable.)  We  use  an  Immersed  Boundary  Fractional  Step  method  with  a 
fast  solver,  developed  by  Colonius  and  Taira  [13j.  The  studies  here  are  for  low  Reynolds 
number  (mostly  at  Re  =  100),  but  the  modeling  techniques  are  equally  applicable  at  higher 
Reynolds  numbers. 
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Figure  1:  Overall  structure  for  coupling  of  flight  dynamics  to  unsteady  aerodynamic  models. 

Lagrangian  Coherent  Structures.  In  order  to  identify  the  coherent  structures  that  arise 
in  these  unsteady  flows,  we  compute  Lagrangian  Coherent  Structures  |24j.  which  precisely 
determine  the  boundaries  of  separation  bubbles  and  other  structures  one  needs  to  model 
in  order  to  accurately  describe  the  unsteady  aerodynamics.  These  methods  originate  from 
the  dynamical  systems  community,  employing  ideas  of  hyperbolic  invariant  manifolds  and 
Finite-Time  Lyapunov  Exponents  (FTLE).  These  techniques  give  a  clear  picture  of  the 
structure  of  these  complex  unsteady  flows,  and  thus  allow  one  to  better  understand  the  flow 
physics. 

Balanced  Proper  Orthogonal  Decomposition.  In  a  related  MURI  project  on  Closed- 
Loop  Flow  Control,  we  have  developed  improved  model  reduction  methods,  based  on  ap¬ 
proximations  of  balanced  truncation  that  are  tractable  even  for  very  large  systems.  This 
method,  called  Balanced  Proper  Orthogonal  Decomposition  (BPOD)  has  been  quite  suc¬ 
cessful  for  producing  low-order  models  suitable  for  control  design  [2],  and  significantly  out¬ 
performs  standard  POD-Galerkin  methods.  In  the  present  effort,  we  employ  this  technique, 
and  develop  it  further,  introducing  a  method  that  does  not  require  adjoint  simulations,  and 
extending  it  to  systems  with  periodic  orbits  (such  as  periodic  vortex  shedding). 

1.2  Summary  of  this  report 

This  report  describes  the  results  of  this  three-year  project,  and  the  technical  results  are 
divided  into  five  sections.  In  Section  [2j  we  provide  some  necessary  background  on  classical 
models  for  unsteady  aerodynamics. 

In  Section  [3j  we  describe  how  we  use  Finite-Time  Lyapunov  Exponents  (FTLE)  to  identify 
when  unsteady  separation  occurs,  and  use  these  quantities  to  gain  an  improved  understand¬ 
ing  of  when  the  classical  models  break  down.  We  demonstrate  that  the  classical  models  fail 
when  there  is  significant  leading-edge  separation,  which  occurs  at  large  Strouhal  numbers  or 
reduced  frequencies.  A  key  result  of  this  section  is  a  method  for  fast  computation  of  FTLE 
fields,  when  one  needs  to  compute  many  fields  at  nearby  times  (as  when  making  a  movie). 

Section  [I]  presents  two  different  modeling  procedures  for  predicting  unsteady  aerodynam¬ 
ics  using  nonlinear  models.  One  approach  is  based  on  Proper  Orthogonal  Decomposition 
(POD)  and  Galerkin  projection,  and  the  other  is  a  simple  “phenomenological”  model  based 
on  the  phenomena  observed  in  the  real  flow.  The  models  include  nonlinear  behavior  such 
as  Hopf  bifurcations  and  limit  cycles,  both  of  which  are  observed  in  the  real  flow.  However, 
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both  modeling  strategies  have  a  serious  limitation  in  that  they  need  to  be  calibrated  for  a 
particular  angle  of  attack,  and  are  not  immediately  suitable  for  control  design.  Neverthe¬ 
less,  they  provide  a  starting  point  for  future  work  in  developing  nonlinear  models  for  these 
unsteady  flows. 

Section  [5]  presents  the  technique  of  model  reduction  that  has  been  most  successful  in  this 
work:  Balanced  Proper  Orthogonal  Decomposition  (Balanced  POD).  This  technique  is  based 
on  linear  models,  and  typically  requires  numerical  simulations  of  both  the  linearized  equa¬ 
tions  and  the  corresponding  adjoint  equations.  The  need  for  adjoint  information  is  a  severe 
limitation  of  this  method,  as  it  prohibits  it  from  being  applied  to  experimental  data.  A  ma¬ 
jor  result  of  this  project  is  a  method  for  developing  reduced-order  models  that  are  identical 
to  Balanced  POD,  without  the  need  for  adjoint  information.  This  procedure  is  presented  in 
Section [5.1 1  Another  extension  is  to  flows  linearized  about  a  periodic  orbit,  such  as  periodic 
vortex  shedding,  and  this  is  described  in  Section  |5.2| 

The  model  reduction  methods  developed  in  Section  [5]  are  then  applied  to  unsteady  aerody¬ 
namics  in  Section  |6[  and  we  view  the  models  in  this  section  as  the  models  most  appropriate 
for  design  of  flight  controllers.  These  models  are  based  on  Wagner’s  indicial  response  [63]. 
and  agree  with  Wagner’s  models  to  an  arbitrarily  high  degree  of  accuracy.  While  Wagner’s 
models  represent  the  lift  as  a  convolution  integral  that  is  cumbersome  to  compute,  and  not 
suitable  for  control  synthesis,  our  procedure  produces  state-space  models  that  may  be  used 
directly  for  control  design.  Furthermore,  our  models  are  formulated  in  a  way  that  directly 
builds  upon  standard  approaches  incorporating  classical  “stability  derivatives”  Cxa,  C^q, 
Clo  (as  in  |57|),  as  shown  in  Figure  [28] 

The  models  obtained  in  Section [6]predict  the  unsteady  response  very  well  for  small-amplitude 
maneuvers.  For  large  amplitude  maneuvers,  agreement  is  poor  when  the  angle  of  attack 
becomes  larger  than  about  20°.  More  work  is  needed,  however,  to  extend  these  models  to 
the  nonlinear  regime. 

2  Background:  classical  models 

When  modeling  the  aerodynamic  forces  acting  on  an  airfoil  in  motion,  it  is  natural  to  start 
with  a  quasi-steady  approximation.  Instead  of  dealing  with  the  full  unsteady  problem,  one 
assumes  that  the  airfoil’s  center  of  mass,  h,  and  angle  of  attack,  a,  motions  are  “gradual” 
enough  for  the  flow  held  to  locally  equilibrate  to  the  motion.  In  this  way,  the  unsteady  terms 
in  the  how  equations  are  set  to  zero  and  the  motion  is  accounted  for  by  translating  h  into 
an  effective  angle  of  attack  and  a  into  an  effective  camber.  Finally,  applying  the  assumption 
of  a  thin  airfoil,  we  obtain  a  quasi-steady  estimate  for  the  lift  coefficient 

Cl  —  27t  h  T 

Lengths  are  nondimensionalized  by  2b  and  time  is  nondimensionalized  by  2b/Uao,  where 
Uoo  is  the  free  stream  velocity,  b  is  the  half-chord  length  and  a  is  the  pitch  axis  location 
with  respect  to  the  1/2-chord  (e.g.,  pitching  about  the  leading  edge  corresponds  to  a  =  — 1, 
whereas  the  trailing  edge  is  a  =  1). 

2.1  Theodorsen’s  frequency  response 

In  1935,  Theodorsen|58j  went  beyond  the  quasi-steady  models  and  solved  for  the  lift  dis¬ 
tribution  around  an  idealized  airfoil  in  purely  sinusoidal  pitch  and  plunge  maneuvers.  His 
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theory  is  soluble  with  analytic  techniques,  relying  on  a  number  of  simplifications,  such  as 
an  inviscid,  incompressible  flow  held  and  infinitesimal  deflections  of  a  hat  plate,  leaving  an 
idealized  planar  wake.  Because  his  theory  was  developed  to  handle  purely  sinusoidal  maneu¬ 
vers,  it  is  represented  in  the  frequency  domain.  Theodorsen’s  model  predicts  the  unsteady 
lift  as 
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where  Theodorsen’s  function  C(k )  is  a  transfer  function  relating  sinusoidal  inputs  of  reduced 
frequency|30|  k  to  their  aerodynamic  response.  The  first  set  of  terms  represent  the  “apparent 
mass”,  or  non-circulatory  terms.  The  second  set  of  terms  are  due  to  circulatory  effects,  and 
are  exactly  the  quasi-steady  forces  multiplied  by  Theodorsen’s  function,  which  accounts  for 
the  change  in  magnitude  and  phase  of  these  terms  with  changes  in  reduced  frequency.  This 
expression  simplifies  considerably  for  an  airfoil  in  pure  pitch  or  pure  plunge 


pure  plunge  (a  =  0) 
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2.2  Wagner’s  indicial  response 

Although  Theodorsen’s  model  is  a  powerful  tool  for  determining  unsteady  lift  coefficients,  it 
is  only  solvable  in  closed  form  for  sinusoidal  forcing.  The  time  domain  method  of  Wagner|63j 
makes  it  possible  to  reconstruct  the  lift  response  to  arbitrary  angle-of-attack  input,  a(t),  by 
superposition  of  the  “indicial”  lift  response  C'£(t)  due  to  a  step  response  in  angle  of  attack, 
a  =  5(t): 

Cl(£)  =  c[{t)a(0)  +  [  c[(t  -  r)d(r)dr  (4) 

Jo 

Wagner  originally  derived  the  indicial  response  analytically,  accounting  for  added  mass  and 
shed-wake  effects  in  a  manor  similar  to  that  of  Theodorsen.  However,  it  is  possible  to  reduce 
the  number  of  simplifying  assumptions  by  obtaining  the  indicial  response  C'£  from  experi¬ 
ment.  Therefore,  the  only  assumption  is  that  of  linearity;  a  more  general  approach  based 
on  functionals  has  been  developed  to  extend  this  theory  for  nonlinear  indicial  response  |59j  ■ 

There  are  a  number  of  interesting  generalizations  to  these  methods,  such  as  the  methods 
of  Sears [62]  and  Kiissner|31|.  which  extended  the  methods  of  Theodorsen  and  Wagner,  re¬ 
spectively,  to  the  problem  of  moving  wind  direction,  rather  than  moving  airfoil.  However, 
this  section  is  meant  to  provide  a  brief  review  of  the  classical  tools  developed  for  modeling 
unsteady  aerodynamic  forces  in  the  1920sM950s.  A  more  complete  treatment  of  the  subject 
can  be  found  in  Leishman[33j. 

Wagner’s  indicial  response  benefits  from  its  time  domain  formulation;  however,  the  superpo¬ 
sition  integral  approach  is  computationally  expensive  and  does  not  fit  nicely  into  a  control 
framework.  In  Section  [6j  we  present  a  systematic  approach  for  obtaining  computationally 
tractable  reduced  order  models  (ROMs)  based  on  the  indicial  response.  Moreover,  these 
models  take  the  form  of  low-dimensional,  state-space  equations,  which  are  ideal  for  the 
application  of  analysis  and  control  techniques. 
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3  Finite-time  Lyapunov  exponents  for  detecting  breakdown  of 
classical  models 


Classical  models  for  unsteady  aerodynamics  often  assume  that  the  flow  remains  attached 
over  a  wing.  Violation  of  this  assumption  is  a  common  reason  for  classical  models  to  break 
down,  and  therefore  detecting  unsteady  separation  is  valuable.  Moreover,  an  understanding 
of  flow  separation  provides  insight  into  the  flow  physics  involved  in  these  phenomena. 


Separation  in  unsteady  flows  is  a  surprisingly  subtle  phenomenon.  In  steady  flows,  the 
criterion  for  separation  is  simple:  at  a  separation  point,  the  shear  stress  vanishes  at  the  walls. 
In  an  unsteady  flow,  this  criterion  is  not  valid,  either  when  applied  to  the  time-averaged 
mean  flow,  or  to  the  instantaneous  unsteady  flow,  and  more  sophisticated  techniques  are 
required  |64| . 


Finite-time  Lyapunov  exponents  (FTLE)  are  a  valuable  tool  for  studying  unsteady  sepa¬ 
ration,  and  identifying  separated  regions.  However,  FTLE  fields  are  reasonably  expensive 
to  compute,  especially  when  one  is  interested  in  an  unsteady  phenomenon,  and  visualizing 
FTLE  fields  at  many  different  times.  In  Section  [3. 1|  we  present  a  technique  for  significantly 
improving  the  speed  of  this  computation.  In  Section  3.2  we  use  this  technique  to  iden¬ 
tify  when  classical  unsteady  models  fail,  and  we  demonstrate  that  this  breakdown  typically 
occurs  when  separation  becomes  significant.  For  further  details  on  these  methods,  see  0®- 


3.1  Fast  computation  of  time-varying  FTLE  fields 

The  theory  and  computation  of  finite-time  Lyapunov  exponents  (FTLE),  also  known  as 
direct  Lyapunov  exponents  (DLE),  is  a  relatively  modern  development  [23,  53].  with  exten¬ 
sions  to  3-dimensional  [20|  |22|  and  n-dimensional  [36]  flows.  FTLE  analysis  has  been  widely 
applied  in  a  number  of  branches  of  fluid  mechanics,  including  fluid  transport  mmm,  bio¬ 
propulsion  [48]  [65],  19] ,  flow  over  airfoils  ]38[  [9]  0 ,  plasmas  @2,  and  geophysical  flows  [34,  35j . 

Because  FTLE  analysis  is  particularly  useful  for  unsteady  flows,  it  is  often  necessary  to 
compute  a  sequence  of  FTLE  fields  in  time  to  visualize  an  unsteady  event.  As  flows  become 
more  complex,  computations  become  increasingly  expensive.  In  particular,  FTLE  calcula¬ 
tions  are  expensive  because  a  large  number  of  particle  trajectories  must  be  integrated  in 
order  to  obtain  a  particle  flow  map,  often  from  stored  velocity  fields.  When  computing  a 
sequence  of  FTLE  fields  in  time,  it  is  possible  to  speed  up  the  computation  considerably 
by  eliminating  redundant  particle  integrations.  One  approach  that  has  been  developed  uses 
adaptive  mesh  refinement  to  reduce  the  number  of  integrations  mmm- 

The  approach  here  is  to  construct  an  approximate  flow  map  by  composing  intermediate 
flow  maps  from  FTLE  held  calculations  at  neighboring  times.  The  first  class  of  how  map 
approximation,  denoted  unidirectional  composition,  constructs  a  how  map  by  composing 
intermediate  how  maps  which  are  all  aligned  in  the  same  time  direction.  The  second  class, 
denoted  bidirectional  composition,  composes  intermediate  how  maps  in  both  positive  and 
negative-time.  The  methods  are  compared  using  analytic  estimates  for  accumulated  error 
and  computation  time  as  well  as  benchmarks  on  a  number  of  example  hows. 

Main  results  In  this  section  we  demonstrate  that  the  unidirectional  method  is  both  fast 
and  accurate,  although  it  requires  significantly  more  memory  than  the  bidirectional  method. 
Orders  of  magnitude  speed-up  may  be  achieved  over  the  standard  method,  and  computa¬ 
tional  improvement  scales  with  the  desired  time  resolution  of  the  FTLE  animation.  The 
bidirectional  method  suffers  from  significant  error  which  is  aligned  with  the  opposite-time 
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Lagrangian  coherent  structures.  To  understand  this  coherent  error,  we  provide  an  error 
analysis  for  both  methods,  and  uncover  an  important  relationship  between  positive-time 
LCS  (pLCS)  and  negative-time  LCS  (nLCS).  In  particular,  in  the  neighborhood  of  a  time- 
dependent  saddle,  particles  near  the  pLCS  flow  into  particles  near  the  nLCS  in  positive 
time. 


3.1.1  Standard  computation  of  FTLE 

Consider  a  time-dependent  velocity  field  u  on  M"  and  a  particle  trajectory  x(t)  which  satisfies 

x  =  u(x(t),t).  (5) 

The  velocity  held,  u,  may  be  an  unsteady  solution  of  the  Navier-Stokes  equation,  although  it 
is  only  assumed  that  u  is  at  least  C°  in  time  and  C 1  in  space.  However,  to  extract  Lagrangian 
coherent  structures  from  the  Hessian  of  the  FTLE  held,  u  must  be  C 2  in  space  [53].  The 
velocity  held  may  be  analytically  defined,  but  is  more  often  obtained  from  experiments  or 
direct  numerical  simulation  which  produce  velocity  held  data  at  discrete  snapshots  over  a 
finite  range  of  time.  A  method  of  computing  hnite-time  Lyapunov  exponents  (FTLE)  on  a 
finite  amount  of  discrete  velocity  held  data  has  been  developed  [23],  !53] . 

Computing  an  FTLE  held  typically  involves  four  steps.  First,  a  grid  of  particles  X$  C  Mn  is 
initialized  over  the  domain  of  interest.  The  particles  are  advected  (i.e. ,  integrated)  with  the 
how  from  initial  time  0  to  final  time  T,  resulting  in  a  time-T  particle  how  map,  defined 
as: 

:  Mn  — »  Mn;  x(0)  i—>  x(0)  +  f  u(x(r),r)dr.  (6) 

Jo 


Next,  the  how  map  Jacobian,  D<F^  is  computed,  usually  by  finite-differencing,  to  obtain  the 
Cauchy-Green  deformation  tensor, 


A  =  (DT^)* 


(7) 


where  *  denotes  transpose.  Finally,  the  largest  eigenvalue,  Amax,  of  this  symmetric  tensor  is 
extracted  and  synthesized  into  an  FTLE  held: 

0  ;  x)  =  j^T  log  yAmaxCAtx)).  (8) 

The  bottleneck  in  this  procedure  is  the  large  number  of  particle  integrations  required  to 
obtain  the  particle  how  map,  Moreover,  if  the  velocity  held  is  time-varying,  it  is 

necessary  to  compute  a  sequence  of  FTLE  helds  in  time  to  visualize  unsteady  events,  as 
shown  schematically  in  Fig.  [2] 

3.1.2  Flow  Map  Approximation 

As  seen  in  Fig.  [2]  the  standard  method  of  computing  a  sequence  of  FTLE  helds  involves 
inefficient  re-integration  of  particles.  The  unidirectional  and  bidirectional  methods  outlined 
below  streamline  the  computation  of  neighboring  FTLE  helds  by  approximating  the  time-T 
how  map,  which  can  be  written  as: 


$(0+T  = 
to 


&N  o  •  •  •  o  d*(2  o  0 

CJV- 1  ci  CO 


(9) 


*lhh+T  ^ 


<k+T  ^ 


®h+T 


■i'.'.  r&B 


;; j>' 


- I-*- 


I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  -  Time 

0  h  2h  3h  ■  ■  ■  T  -  -  - 
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Essential 
— >•  Exact  Flow 


Figure  2:  The  standard  method  for  computing  FTLE.  Exact  flow  maps  for  A;  E 

{0,1, 2, 3}  are  shown  (solid  black  arrow).  Essential  (blue)  and  redundant  (red)  particle 
integrations  are  outlined  in  dashed  ovals. 


where  f/v  =  to  +  T. 

Because  the  flow  maps  are  obtained  numerically  on  a  discrete  grid  of  points,  Xq  C  Mn,  it 
is  necessary  to  interpolate  the  map  at  points  x  ^  Xq.  Consider  a  flow  map  T  :  Mn  — ►  Mn, 
defined  on  Mn,  and  the  same  flow  map  restricted  to  Xq.  d>|x0  :  Xq  — ►  Mn.  The  interpolation 
operator  X  takes  the  discrete  map  T|x0  and  returns  the  map,  T,  defined  on  all  of  Mn: 

X  :  $\Xo  ^  T  (10) 

Using  the  shorthand  XT  =  X (T|x0)>  we  obtain  the  approximate  flow  map: 

ig+T(Vo)  =  !*;»_,  O  •  •  •  O  l<S>t  O  4 j;(X„) 

»  *£+T(*o) 

The  bidirectional  method  approximates  the  time-T  flow  map  by  first  integrating  back¬ 

ward  to  a  reference  time,  t  =  0,  then  interpolating  forward  through  a  previously  computed 
time-T  map,  T^,  and  finally  integrating  forward  to  time  to  +  X.  The  unidirectional  method 
approximates  the  time-T  flow  map  using  a  number  of  smaller  time  flow  maps,  which 

all  have  the  same  time  direction.  Additionally,  the  chain  rule  may  be  applied  to  each  of  the 
methods,  resulting  in  an  approximation  to  the  flow  map  Jacobian,  DT^+T. 

3.1.3  Bidirectional  Composition 

Bidirectional  approximation  eliminates  redundancy  from  neighboring  FTLE  field  computa¬ 
tions  by  using  the  information  from  a  known  flow  map  at  a  given  time,  Tq  ,  to  calculate  an 
approximation  to  the  flow  map  at  future  times,  First,  Xq  is  integrated  backward 

from  to  to  the  reference  time  0.  The  distorted  grid  T®0(Ao)  is  then  flowed  forward  through 
the  interpolated  map,  XTq  ,  and  finally  integrated  forward  an  amount  to  to  the  desired  time 
to  +  T,  as  in  Fig.  [3] 

T|;°+T  =  ^0+ToX$[o$“i0.  (12) 

The  flow  is  stored  as  a  reference  solution  to  compute  an  approximation  to  the  flow  map 
at  later  times  T^+t  fs  T^+t  by 

T^+t  =  h+T  o  XT^  keZ  (13) 
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Figure  3:  Schematic  for  bidirectional  method  (a).  Given  a  known  flow  map  (solid  black 
arrow),  it  is  possible  to  approximate  the  flow  map  at  later  times  (dashed  black 

arrow)  by  integrating  backward  in  time  to  t  =  0  (red  arrow),  flowing  forward  through  the 
interpolated  map  ZTq  which  was  already  computed  (blue  double  arrow),  and  integrating 
trajectories  forward  to  the  correct  final  time  (green  arrow). 


Instead  of  using  <Fq  as  the  reference  solution  for  every  future  time,  it  is  convenient  to  use 
the  new  approximate  flow  map  as  the  reference  solution  for  the  next  iteration: 


(14) 


This  method  may  be  continued,  using  to  approximate 


; Ak+i)h+T . 


$ 


(k+l)h.+T  _  ^(fc+ljh+T  ■jx.k} i+T  x.kh 

(k+l)h  ~Vkh+T  °-LVkh  °^(k+l)h- 


(15) 


Errors  will  compound  more  quickly  since  we  are  using  approximate  flow  maps  as  the  reference 
solutions  for  later  approximations,  as  seen  in  Fig.  [f] 


H — I — I — I — I — I — I — I — 1 — h 


HI - 1 - 1 - 1 - H ►  Time 


Known  Flow  Map 


Figure  4:  Schematic  for  bidirectional  method  (b).  As  in  Fig.  [3j  a  known  flow  map  (solid 
black  arrow)  is  used  to  approximate  the  flow  map  at  a  later  time  <f>^+T  (dashed  black 
arrow).  The  approximate  flow  map  is  used  as  the  known  map  for  the  next  step  (dashed 
black  arrow). 


3.1.4  Unidirectional  Composition 

The  basis  of  the  unidirectional  method  is  to  eliminate  redundant  particle  integrations  by  only 
integrating  particle  trajectories  through  a  given  velocity  field  a  single  time.  If  a  sequence  of 
FTLE  snapshots  is  desired  at  a  time  spacing  of  h,  for  example  as  frames  in  an  animation, 
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then  it  is  convenient  to  break  up  the  time-T  flow  map  into  smaller  time-h  flow  maps,  where 
T  =  kh: 


*Qh  =  l^_x)ho-..ol^hho^  (16) 

This  method  is  called  unidirectional  because  particle  flow  maps  of  the  same  time  direction 
are  used,  as  opposed  to  the  bidirectional  method  which  composes  both  positive-time  and 
negative-time  flow  maps. 


0  h  2h  3h  ■  ■  ■  T  T+h  ■  ■  ■ 


Desired  Flow  Map 

Figure  5:  Schematic  for  unidirectional  method.  Time-/i  flow  maps  (short  blue  arrows)  are 
stored  and  composed  to  approximate  the  time-T  flow  map  (long  black  arrow).  The  next 
flow  map  only  requires  integrating  one  new  time-h  flow  map  (green  arrow). 

The  simplest  approach  is  to  compute  a  number  of  tirne-h  flow  maps  and  store  them  in 
memory.  Then,  to  construct  an  approximate  it  remains  only  to  compose  the  sequence 

of  interpolated  time-/i  flow  maps.  The  next  iteration  involves  integrating  one  more  tirne-h 
flow  map  and  composing  the  next  sequence,  as  in  Fig.  [5] 

To  further  improve  efficiency  by  reducing  the  total  number  of  flow  map  compositions,  it  is 
possible  to  construct  a  multi-tiered  hierarchy  of  flow  maps  for  reuse  in  neighboring  flow  map 
constructions.  Given  enough  memory,  it  is  possible  to  reduce  the  number  of  interpolated 
compositions  by  increasing  the  number  of  tiers  of  flow  maps,  each  tier  being  constructed  as 
the  composition  of  two  of  the  flow  maps  in  the  next  tier  lower,  as  in  Fig.  [6] 


->!- 


I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  »  Time 

0  h  2h  3h  ■  ■  ■  T  h+T  ■  ■  ■ 


Desired  Flow  Map 


Figure  6:  Schematic  for  unidirectional  method  with  multiple  tiers.  The  bottom  tier  of  time- 
h  flow  maps  is  computed  as  in  Fig.  [5}  Pairs  are  composed  to  form  the  next  tier  of  time-2/i 
flow  maps,  and  so  on.  This  method  requires  more  storage,  but  fewer  total  compositions 
when  computing  a  series  of  FTLE  fields  for  an  animation. 


3.1.5  Chain  Rule  of  Compositions 

As  seen  in  Eq.  (|t|)  ,  once  the  flow  map  is  obtained,  it  is  necessary  to  compute  the  flow 

map  Jacobian  in  order  to  extract  the  FTLE.  Applying  the  chain  rule  to  Eq.  it  is  possible 
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Problem 

Resolution 

T/h 

Frames 

Method 

Mem.  (GB) 

Speed-up 

Accurate 

Double  Gyre 

1024  x  512 

15 

30 

Standard 

.05 

1 

Yes 

Unidirectional 

.36 

10 

Yes 

Bidirectional 

.14 

6.2 

No 

Pitching  plate 

1024  x  512 

15 

30 

Standard 

.48 

1 

Yes 

Unidirectional 

.70 

8.2 

Yes 

Bidirectional 

.50 

5.4 

No 

Pitching  plate 

600  x  300 

150 

192 

Standard 

.48 

1 

Yes 

Unidirectional 

1.8 

67 

Yes 

Bidirectional 

.48 

54 

No 

ABC  flow 

1283 

20 

40 

Standard 

.48 

1 

Yes 

Unidirectional 

2.6 

6.8 

Yes 

Bidirectional 

.73 

7.3 

No 

Table  1:  Comparison  of  methods  on  various  examples  fluid  flows.  The  unidirectional  method 
both  fast  and  accurate,  but  requires  more  memory  than  the  other  methods,  providing  one 
or  two  orders  of  magnitude  computational  improvement  over  the  standard  method. 


to  express  the  flow  map  Jacobian  as  a  product  of  the  Jacobians  of  intermediate  flow  maps: 


D(^)(x)  =  D(T 


to 


tN-l 


= 

tN-l 


$ 


tN-l  i 


to  (x)) 


(x) 

x-xD$J(x) 


Applied  to  the  bidirectional  methods,  this  yields: 


$h+T  =$h+T 


0  $o  °  K 


=►  D^^T(x)  =D<^+t  K  °  *°)  (x)  x 

xD^K)  (x)oDlJ(x), 

and  applied  to  the  unidirectional  methods,  this  yields: 


<  =®T-M 
D^(x)  =D 


o  $lh  °  $0 


-JT—h 

■0 


•  •  x  D $lh 


(x))  X- 
T0h(x))  x  DTg(x). 


(17) 


(18) 


(19) 


3.1.6  Comparison  of  Methods 

Each  method  from  Section[3.1.2  is  implemented  and  tested  on  three  example  problems:  the 
periodic  double  gyre,  2D  flow  over  a  pitching  flat  plate  at  Reynolds  number  100,  and  3D 
unsteady  ABC  flow.  These  examples  are  chosen  because  they  cover  a  range  of  features 
including  2D  and  3D  vector  fields,  which  are  either  defined  analytically  or  obtained  from 
data  files  from  DNS  on  either  open,  closed,  or  periodic  domains.  Each  example  problem 
is  discussed  more  in  Appendix  B  of  [5],  including  details  such  as  how  the  velocity  field  is 
defined,  and  on  what  domain.  In  the  pitching  plate  example,  velocity  field  snapshots  are  all 
loaded  up-front  before  applying  the  methods. 
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(Cl)  Exact  (C2)  Unidirectional  (C3)  Bidirectional  (C4)  Exact  (opposite-time) 


(Rl)  Double  Gyre 


(R2)  Pitching  Plate 


(R3)  Pitching  Plate 


(R4)  Unsteady  ABC 


Figure  7:  Graphical  comparison  of  each  method  on  four  examples:  (top  row)  positive-time 
FTLE  of  double  gyre,  (second  row)  positive-time  FTLE  of  2D  pitching  plate,  (third  row) 
negative-time  FTLE  of  2D  pitching  plate,  (bottom  row)  negative-time  FTLE  of  3D  ABC 
flow.  Each  figure  shows  the  FTLE  held  after  a  number  of  iterations  of  the  given  method.  The 
column  of  FTLE  fields  calculated  using  unidirectional  composition  agree  well  with  the  exact 
FTLE  fields  computed  using  the  standard  method.  The  column  of  FTLE  fields  calculated 
using  bidirectional  composition  all  have  significant  error  which  is  aligned  with  the  opposite¬ 
time  coherent  structures.  The  opposite-time  FTLE  fields  are  shown  in  the  rightmost  column 
for  comparison  with  the  bidirectional  method.  FTLE  fields  computed  for  positive-time  how 
maps  are  blue  and  those  computed  for  negative-time  how  maps  are  red. 
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Figure  8:  Convergence  tests  for  FTLE  field  error  vs.  integration  time-step  and  grid  spacing 
on  double-gyre. 


Table  [l]  summarizes  the  results  comparing  each  method  on  the  three  example  fluid  flows. 
In  each  comparison,  the  standard,  unidirectional  and  bidirectional  methods  are  used  to 
compute  a  sequence  of  FTLE  fields  which  are  frames  in  an  unsteady  animation.  The  flow 
map  duration  used  to  compute  an  FTLE  field  is  T,  and  the  time-spacing  between  neighboring 
FTLE  fields  is  h,  so  the  number  of  animation  frames  per  flow  map  duration  is  T/h.  As 
demonstrated  in  Section  3.1.8|  this  is  an  upper  bound  on  the  speed-up  of  the  unidirectional 
method. 


In  each  comparison,  the  unidirectional  method  is  accurate  and  offers  the  greatest  speed-up 
over  the  standard  method.  However,  it  also  requires  more  memory  than  any  other  method. 
The  bidirectional  method  is  fast  and  uses  less  memory  than  the  unidirectional  method,  but 
is  prone  to  large  errors  in  the  approximate  flow  map  and  does  not  accurately  reproduce  the 
FTLE  field. 

Contour  plots  of  the  FTLE  fields  computed  after  a  number  of  iterations  of  each  method 
are  shown  in  Fig.  [7J  The  FTLE  fields  computed  with  the  unidirectional  method  agree  with 
those  computed  using  the  standard  method,  as  seen  by  comparing  the  first  and  second 
column  of  Fig.  [7]  FTLE  fields  computed  using  the  bidirectional  method,  shown  in  the  third 
column,  have  large  errors.  It  is  interesting  to  note  that  these  errors  are  aligned  with  coherent 
structures  found  in  the  opposite-time  FTLE  field,  shown  in  the  fourth  column.  An  analysis 
of  this  coherent  error  is  provided  in  Section  V  of  0. 

3.1. 7  Example  -  Double  Gyre 

Figure  [8]  shows  the  L2  and  La 0  error  of  the  forward-time  FTLE  field  for  the  double  gyre 
computed  using  the  standard  method  with  T  =  16,  as  time-step  At  and  grid  spacing  Ax 
are  varied.  At  a  given  grid  spacing,  a  reference  FTLE  field  is  computed  using  a  sufficiently 
small  time-step,  At  =  10”4,  so  that  the  FTLE  field  may  be  considered  exact.  For  small 
enough  time-step  At  ~  .001,  the  FTLE  field  error  converges.  All  integrations  are  performed 
using  a  fixed  time-step,  fourth  order  Ruuge-Kutta  scheme. 
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Figure  9:  Computational  time  vs.  Iteration  (top)  and  L2  error  vs.  Iteration  (bottom)  for 
FTLE  fields  of  the  double  gyre  with  resolution  1024  x  512. 


The  flow  map  approximation  methods  are  only  faster  than  the  standard  method  when  used 
to  compute  a  sequence  of  FTLE  fields  in  time,  as  in  the  construction  of  frames  for  a  movie. 
Figure  [9]  compares  computation  time  and  L2  error  vs.  frame  number  (iteration  jf)  for  a 
sequence  of  FTLE  fields  of  the  double  gyre,  computed  using  the  standard,  unidirectional, 
and  bidirectional  methods.  Each  iteration  produces  an  FTLE  field  which  is  a  single  frame 
in  an  animation  of  the  unsteady  FTLE  field.  In  this  example,  the  flow  map  duration  is 
T  =  16,  the  time  spacing  between  each  FTLE  field  is  h  =  1,  and  the  time-step  of  integration 
is  At  =  .01.  The  multi-tier  unidirectional  method  uses  four  tiers. 

The  first  FTLE  field  takes  roughly  the  same  time  to  compute  using  each  of  the  meth¬ 
ods.  However,  for  subsequent  iterations,  the  unidirectional  and  bidirectional  methods  are 
significantly  faster.  The  computation  time  of  bidirectional  method  (a)  increases  with  the 
number  of  iterations,  k,  because  integrating  back  from  t  =  kh  to  the  reference  time  t  =  0 
becomes  more  costly  as  k  increases,  as  seen  in  Fig.  [d]  After  T/2/i  =  8  iterations  of  bidi¬ 
rectional  method  (a),  it  is  advantageous  to  compute  a  new  reference  flow  map  using  the 
standard  method.  This  explains  the  breaks  in  the  solid  red  curve  in  part  (b)  of  Fig.  [9]  as  the 
bidirectional  method  is  exact  at  these  iterations.  Bidirectional  method  (b)  overcomes  this 
increasing  cost  vs.  iteration  by  using  the  flow  map  from  the  current  iteration  as  the  reference 
flow  map  at  the  next  iteration.  However,  using  an  approximate  flow  map  to  compute  the 
next  approximation  causes  bidirectional  method  (b)  to  accumulate  error  more  quickly  than 
method  (a).  The  unidirectional  method  is  both  the  fastest  and  most  accurate  method  in 
this  comparison. 

3.1.8  Computational  Resources 

Again,  consider  a  sequence  of  time-T  flow  maps  spaced  h  apart,  as  might  be  required  for  an 
unsteady  visualization.  When  there  are  many  integration  time- steps  of  size  At  between  each 
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neighboring  flow  map,  i.e.  At  <C  h,  then  the  added  cost  of  flow  map  composition  becomes 
relatively  small  compared  with  the  cost  of  integrating  a  time- It  flow  map. 

The  standard  method  involves  ( T/h )  X  {h/ At)  integration  steps  for  each  new  FTLE  field, 
whereas  the  unidirectional  method  only  requires  h/ At  integration  steps,  and  the  bidirec¬ 
tional  method  requires  2h/At  integration  steps.  Assuming  At  <C  h.  the  speed-up  of  the 
unidirectional  method  over  the  standard  method  will  increase  as  the  number  of  frames  in 
the  animation  per  flow  map  duration,  T.  In  other  words,  as  A  t/h  — ►  0,  the  computation  of 
using  the  unidirectional  method  is  T/h  times  faster  than  using  the  standard  method, 
and  twice  as  fast  as  the  bidirectional  method. 

In  the  examples  above,  all  intermediate  flow  maps  were  stored  in  memory  until  no  longer 
useful  for  future  computations.  Regardless  of  any  parameters  of  the  FTLE  field  animation, 
the  standard  and  bidirectional  methods  must  store  a  fixed  number  of  flow  maps.  The 
standard  method  stores  the  single  flow  map  while  the  bidirectional  method  stores 

three  maps:  <Lj0,  and  <&t/'i+T .  The  unidirectional  method,  however,  stores  of  every 
intermediate  time-/i  flow  map  &f/_  of  which  there  are  T/h.  Therefore,  the  memory 
requirement  of  the  unidirectional  method  scales  linearly  with  the  upper-bound  on  its  speed¬ 
up,  T/h. 

The  memory  usage  of  the  unidirectional  method  scales  with  the  dimension  of  the  flow  D, 
the  spatial  resolution  R ,  and  the  possible  computational  speed  up  of  the  method  S,  given 
by  T/h : 


Memory  (GB)  ~  S  x  D  x  RD 


8  B/double 
10243  B/GB 


x  flx  Rd 


(20) 

(21) 


For  example,  a  series  of  two  dimensional,  high-definition  (1920  X  1080  resolution)  FTLE 
fields  may  be  computed  using  the  unidirectional  method  with  up  to  100  X  speed  up  using 
approximately  3.1  GB  of  RAM.  A  three  dimensional  FTLE  field  with  resolution  512x256x64 
may  be  computed  with  up  to  100 x  speed  up  with  approximately  19  GB  of  RAM. 

In  the  pitching  plate  example,  velocity  fields  are  obtained  from  data  files  which  are  the  output 
of  a  direct  numerical  simulation.  Because  loading  velocity  fields  which  are  stored  on  disk  is 
slow,  it  is  important  to  minimize  the  number  of  file  loads.  In  the  pitching  plate  example,  all 
of  the  velocity  fields  are  loaded  up-front  and  stored  in  memory  throughout  the  computation. 
However,  velocity  fields  are  often  too  large  to  store  them  all  in  memory,  for  example  in  large 
2D  or  3D  simulations.  After  the  unidirectional  method  is  initialized,  subsequent  iterations 
of  the  method  only  require  loading  velocity  fields  involved  in  the  computation  of  a  single, 
time-/i  flow  map.  The  standard  method,  however,  must  load  velocity  fields  relevant  to  the 
entire  time-T  flow  map,  requiring  T /h  times  as  many  file  loads  as  the  unidirectional  method. 

3.2  Breakdown  of  classical  unsteady  models 

Here,  we  examine  the  effectiveness  of  the  classical  models  described  in  Section [2]  for  an  airfoil 
undergoing  sinusoidal  pitching  or  plunging.  The  simulations  here  are  for  a  flat  plate  at  low 
Reynolds  number  (Re  =  100),  and  are  calculated  using  the  immersed  boundary  method 
of  [15],  For  details  of  the  simulation,  see  [7]. 

For  this  low-Reynolds-number  flat  plate,  the  steady-state  lift  is  shown  as  a  function  of  angle 
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Figure  10:  Lift  coefficient  vs.  angle  of  attack  for  fixed  plate  at  Re  =  100.  A  supercritical 
Hopf  bifurcation  occurs  at  a.  ~  28°.  The  dashed  line  represents  average  post-stall  lift  and 
the  dotted  lines  represent  minimum  and  maximum  post-stall  lift. 


of  attack  in  Figure  10  At  ac  ~  28°,  a  Hopf  bifurcation  occurs[2],  and  for  a  >  ctc,  the  flow 


is  unsteady,  corresponding  to  periodic  vortex  shedding  from  the  leading  and  trailing  edges. 


3.2.1  Theodorsen’s  model  —  sinusoidally  plunging  airfoil 

Using  equations  ([I])  and  (|2])  it  is  possible  to  compare  the  thin  airfoil  theory  and  Theodorsen’s 
model  with  the  simulated  response  of  a  flat  plate  to  sinusoidal  plunging  in  Reynolds  number 
100  flow.  It  is  also  compared  with  an  effective  angle  of  attack  approximation  using  a  look  up 
table  for  the  lift  coefficient  at  the  static  angle  of  attack  aefr(t)  =  tan~1(—h(t)/UOQ)  for  each 
point  in  the  maneuver;  this  approximation  may  be  classified  as  quasi-steady.  The  plunging 
motion  is  specified  by  the  center  of  mass  motion  h(t)  =  —  Asm(u)t). 


It  has  been  shown  [9j  that  for  reasonable  Strouhal  numbers  and  reduced  frequencies  k  = 
irfc/Uoo  less  than  0.5,  the  effective  angle  of  attack  approximation  agrees  well  with  DNS. 
However,  for  the  same  range  of  Strouhal  numbers  St  £  {.032,  .064,  .128},  we  find  that 


Theodorsen’s  theory  agrees  with  DNS  up  to  reduced  frequencies  of  2.0,  shown  in  Figure  11 


This  is  consistent,  since  at  larger  reduced  frequencies,  quasi-steady  assumptions  break  down 
and  it  becomes  important  to  consider  flow  acceleration  terms. 


For  larger  Strouhal  numbers,  St 


even  for  small  reduced  frequencies,  shown  in  Figure  12 


£  {.256,  .512},  Theodorsen’s  model  disagrees  with  DNS 

This  is  particularly  interesting, 


because  these  large  Strouhal  numbers  correspond  to  maximum  effective  angles  of  attack 
which  are  larger  than  the  critical  stall  angle  shown  in  Figure  [TO]  At  these  Strouhal  num¬ 
bers  the  effective  angle  of  attack  approximation  plateaus  due  to  aeff  >  ac ■  In  addition  to 
disagreeing  with  DNS  in  magnitude  and  phase,  Theodorsen’s  model  does  not  describe  the 
higher-frequency  components  in  the  response. 


The  table  below  shows  the  maximum  effective  angle  of  attack  associated  with  each  Strouhal 
number  used  above.  Notice  that  for  Strouhal  numbers  .256  and  .512,  the  maximum  effective 
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St  =  .032,  Reduced  Frequency  k  =  1.0  St  =  .032,  Reduced  Frequency  k  =  4.0 


Heaving  w/  A=.05,  ca=2.0 


Heaving  w/  A=.01 ,  m=8.0 


Figure  11:  Each  curve  is  a  plot  of  lift  coefficient  vs.  time  for  a  sinusoidally  plunging  flat 
plate  at  Re  =  100  and  Strouhal  number  .032.  The  blue  curve  is  the  Cl  from  DNS,  the 
red  curve  is  computed  using  thin  airfoil  theory,  Eq.  ([!]),  the  black  curve  uses  an  effective 
angle  of  attack  approximation,  and  the  green  curve  is  Theodorsen’s  prediction,  Eq.  (J2|. 
Theodorsen’s  model  agrees  well  with  DNS  for  reduced  frequencies  k  <  2.0  as  long  as  the 
Strouhal  number  is  small  enough  that  the  maximum  effective  angle  of  attack  is  less  than  the 
stall  angle. 


St  =  .256,  Reduced  Frequency  k  =  .25  St  =  .512,  Reduced  Frequency  k  =  .5 


Heaving  w/  A=1 .6,  ca=0.5 


Heaving  w/  A=1 .6,  co=1 .0 


Figure  12:  Each  curve  is  a  plot  of  lift  coefficient  vs.  time  for  a  sinusoidally  plunging  flat 
plate  at  Re  =  100  and  Strouhal  numbers  .256  and  .512.  The  blue  curve  is  the  Cl  from 
DNS,  the  red  curve  is  computed  using  thin  airfoil  theory,  Eq.  ([Tj) ,  the  black  curve  uses  an 
effective  angle  of  attack  approximation,  and  the  green  curve  is  Theodorsen’s  prediction,  Eq. 
©•  For  large  Strouhal  numbers,  Theodorsen’s  model  does  not  agree  well  with  DNS  even  at 
low  reduced  frequency.  Also,  there  are  higher  frequency  components  which  are  not  captured 
by  the  models. 


18 


Strouhal  Number  St  =  .274 


Strouhal  Number  St  =  .137 


Pitching  w/a=20,  f=0.4,  leading  edge 


Pitching  w/  a=20,  f=0.4,  mid  point 


Time 


Time 


Figure  13:  Each  curve  is  a  plot  of  lift  coefficient  vs.  time  for  a  sinusoidally  pitching  flat 
plate  at  Re  =  100  and  reduced  frequency  k  =  1.26.  The  blue  curve  is  the  Cl  from  DNS,  and 
the  green  curve  is  Theodorsen’s  prediction,  Eq.  ([3]).  Theodorsen’s  model  matches  DNS  for 
a  flat  plate  pitching  to  an  amplitude  of  20°  about  the  leading  edge  and  not  the  half-chord, 
even  though  the  Strouhal  number  is  larger  for  pitching  about  the  leading  edge. 


angle  of  attack  is  larger  than  the  critical  stall  angle  ac  ~  28°. 


Strouhal  number  (St) 

.032 

.064 

.128 

.256 

.512 

Max  effective  aoa  (a™x) 

5.74° 

11.37° 

21.91° 

38.81° 

58.13° 

3.2.2  Theodorsen’s  model  —  sinusoidally  pitching  airfoil 

Similar  to  the  case  of  a  sinusoidally  plunging  airfoil,  with  equation  ([3])  it  is  possible  to 
compare  Theodorsen’s  model  with  DNS  for  a  sinusoidally  pitching  flat  plate.  Three  pitching 
amplitudes  a  £  {20°,  27.1°,  43.2°}  are  examined,  and  for  each  amplitude  the  pitching  point 
is  varied  along  the  chord  from  the  leading-edge  to  the  trailing-edge  at  every  quarter-chord 
in  between.  In  each  case,  the  reduced  frequency  is  k  =  1.26.  Interestingly,  the  agreement 
between  Theodorsen’s  model  and  DNS  does  not  depend  so  much  on  reduced  frequency  and 
Strouhal  number  as  in  the  case  of  sinusoidal  plunging,  but  instead  depends  much  more  on 
raw  pitching  amplitudes  and  pitching  point.  In  Figure [l3| we  see  that  for  the  same  angle  of 
attack  excursion  Theodorsen’s  model  agrees  better  with  DNS  when  the  plate  pitches  about 
the  leading  edge,  even  though  the  Strouhal  number  is  larger. 

This  is  consistent  with  the  fact  that  Theodorsen’s  theory  depends  on  the  assumption  of 
attached  flow  over  the  wing  surface,  and  pitching  about  the  mid-chord  promotes  leading- 
edge  separation  much  more  than  pitching  about  the  leading  edge.  This  effect  is  shown  in 
Figure  |TT[  where  the  effect  of  leading  edge  separation  is  exaggerated  due  to  the  large  pitch 
amplitude  amax  =  27.1°. 


When  the  pitch  amplitude  is  larger  than  the  stall  angle  ac  ~  28°,  Theodorsen’s  model  does 


not  agree  with  DNS  even  when  the  plate  is  pitched  about  the  leading  edge;  see  Figure  15 


This  is  interesting  because  there  is  a  similar  observation  in  the  case  of  sinusoidal  plung- 
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Strouhal  Number  St  =  .182 


Time 


Figure  14:  Leading  edge  separation  is  visualized  using  FTLE  fields  for  a  plate  pitching  to 
an  amplitude  of  27.1°  about  the  mid  chord. 


ing,  where  Theodorsen’s  model  begins  to  disagree  at  Strouhal  numbers  which  correspond 
to  effective  angles  of  attack  larger  than  ac.  This  should  not  be  surprising,  however,  since 
Theodorsen’s  theory  is  developed  for  infinitesimal  oscillations  where  the  flow  is  never  sepa¬ 
rated  and  the  wake  is  assumed  to  be  planar.  Under  these  assumptions,  Theodorsen’s  model 
cannot  possibly  take  high  angle  of  attack  vortex  shedding  effects  into  account. 


In  addition  to  the  lack  of  agreement  between  Theodorsen’s  model  and  DNS  for  large  pitching 
amplitudes,  there  are  Cl  variations  at  twice  the  pitching  frequency,  as  there  was  in  the  case 
of  sinusoidal  plunging  for  Strouhal  numbers  St  G  {.256,  .512}.  Using  FTLE  to  visualize  the 
flow  structures,  it  is  possible  to  see  not  only  leading  edge  separation,  but  also  a  distortion 
of  the  FTLE  near  the  plate  near  the  mid  chord,  as  shown  in  Figure  16  It  is  likely  that 


natural  vortex  shedding  at  this  high  angle  of  attack  is  interacting  with  the  separation  due 
to  the  airfoils  pitching  motion. 


3.2.3  Indicial  response 


Compared  with  Theodorsen’s  method  of  predicting  response  forces,  indicial  response,  equa¬ 
tion  Q  is  an  empirical  method  which  relies  on  knowing  only  the  response  in  lift  to  a  small 
step  in  angle  of  attack.  For  the  simulations  below,  the  step  in  a  was  approximated  by  a 
steep  sigmoidal  step  of  1°.  The  indicial  response  roughly  predicts  the  initial  peak  observed 
in  DNS  for  fast  pitch-up  maneuvers  of  moderate  amplitude,  a  =  8°  and  a  =  16°,  shown  in 
Figure  17  However,  superposition  of  a  number  of  small  steps  fails  to  reproduce  transient 


oscillations  as  the  initial  peak  dies  off. 


For  larger  angle  of  attack  pitch-up  maneuvers,  say  a  =  32°,  the  indicial  response  predicts 
the  rough  form  of  transient  lift,  but  doesn’t  capture  the  jagged  peak  which  is  observed  in 
DNS.  Because  the  method  of  indicial  response  involves  staggered  superposition  of  a  number 
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Strouhal  Number  St  =  .548 


Strouhal  Number  St  =  .274 


Pitching  w/a=43.2,  f=0.4,  leading  edge 


Pitching  w/a=43.2,  f=0.4,  mid  chord 


Figure  15:  Each  curve  is  a  plot  of  lift  coefficient  vs.  time  for  a  sinusoidally  pitching  flat  plate 
at  Re  =  100  and  reduced  frequency  k  =  1.26.  The  blue  curve  is  the  Cl  from  DNS,  and  the 
green  curve  is  Theodorsen’s  prediction,  Eq.  (J3j).  For  pitching  amplitude  amax  =  43.2°  >  ac, 
Theodorsen’s  model  doesn’t  agree  with  DNS,  despite  the  pitching  point. 


Time 


Figure  16:  Leading  edge  separation  is  visualized  using  FTLE  fields  for  a  plate  pitching  to 
an  amplitude  of  43.2°  about  the  mid  chord. 
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Pitch  Upw/a=8  ,  duration=.16  Pitch  Up  w/a=16  ,  duration=.16 


Figure  17:  (top)  Comparison  of  lift  coefficient,  Cl,  between  DNS  and  Indicial  response  for 
a  flat  plate  in  pitch-up  maneuver  at  Re  =  300.  (bottom)  Angle  of  attack  vs.  time,  (left) 
Pitch-up  of  8°  with  duration  of  16  time  steps,  (right)  Pitch-up  of  16°  with  duration  16  time 
steps 


of  small  steps  to  reconstruct  a  large  step,  it  is  not  possible  for  this  method  to  predict  the 
periodic  vortex  shedding  which  is  characteristic  of  high  angles  of  attack.  Therefore,  even  if 
it  is  able  to  predict  transient  lifts,  it  is  not  useful  for  steady  state  prediction  at  large  angle 


of  attack,  as  shown  in  Figure  18 


3.2.4  Conclusions 


In  this  section,  we  have  investigated  the  unsteady  aerodynamic  forces  on  low-Reynolds 
number  wings  at  high  angle  of  attack  and  in  pitch  and  plunge  maneuvers  using  2D  direct 
numerical  simulations.  The  classical  theories  of  Theodorsen  and  Wagner  have  been  compared 
with  DNS  for  a  number  of  pitch  and  plunge  maneuvers  of  varying  Strouhal  number,  reduced 
frequency,  pitch  amplitude  and  center.  In  addition  to  determining  when  these  theories 
break  down,  the  flow  field  is  investigated  using  FTLE  to  visualize  relevant  flow  structures 
to  determine  how  the  theories  break  down,  indicating  possible  improvements  to  the  models. 

Comparison  of  Theodorsen’s  model  for  the  lift  of  a  sinusoidally  plunging  flat  plate  with 
forces  from  DNS  showed  agreement  for  moderate  reduced  frequencies  k  <  2.0  for  a  range  of 
Strouhal  numbers  for  which  the  maximum  effective  angle  of  attack  is  smaller  than  the  critical 
stall  angle.  For  the  case  of  a  sinusoidally  pitching  plate,  agreement  between  Theodorsen’s 
model  and  DNS  was  less  dependent  on  Strouhal  number  than  the  position  of  the  pitch  axis 
along  the  chord.  Pitching  about  the  mid-chord,  while  resulting  in  a  smaller  Strouhal  number 
than  pitching  about  the  leading  edge,  promotes  leading  edge  separation  and  dynamic  stall 
effects  which  are  not  captured  by  Theodorsen’s  model.  However,  Theodorsen’s  model  does 
not  agree  with  DNS  for  any  pitch  point  if  the  angle  of  attack  excursion  is  large  enough  to 
cause  periodic  vortex  shedding.  This  is  an  important  relationship  between  the  theory  for 
pitching  and  plunging  airfoils;  in  particular,  the  theory  breaks  down  in  both  cases  when  the 
angle  of  attack  (resp.  effective  angle  of  attack)  excursion  exceeds  the  critical  stall  angle. 
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Pitch  Up  w/  a=32  ,  duration=.64 


Figure  18:  (left)  Comparison  of  lift  coefficient,  Cl,  between  DNS  and  Indicial  response  for 
a  flat  plate  in  pitch-up  at  Re  =  300  and  sinusoidal  pitch  maneuver  at  Re  =  100.  (right) 
FTLE  field  visualization  of  the  periodic  laminar  vortex  shedding  which  takes  place  after  the 
transients  die  down. 


This  observation  is  supported  by  the  method  of  indicial  response,  where  agreement  between 
model  and  data  begins  to  break  down  for  large  pitch-up  maneuvers.  The  inability  to  capture 
unsteady  effects  due  to  high  angle  of  attack  is  a  fundamental  limitation  of  both  methods. 


4  Nonlinear  models  for  unsteady  flows  at  fixed  angle  of  attack 


A  first  step  towards  developing  models  for  unsteady  aerodynamics  is  to  develop  models  valid 
at  a  fixed  angle  of  attack.  In  this  section,  we  summarize  our  results  using  two  approaches.  In 
the  first  approach,  we  develop  phenomenological  models,  which  capture  the  essential  features 
observed,  namely  a  Hopf  bifurcation  at  a  critical  angle  of  attack,  as  shown  in  Figure  10  In 
this  model,  the  angle  of  attack  a  appears  explicitly,  so  it  is  straightforward  to  see  how  it 
would  generalize  it  to  other  angles  of  attack.  However,  in  principle,  such  a  generalization  is 
difficult,  as  the  model  contains  three  empirical  parameters  which  would  need  to  be  adjusted. 


In  the  second  approach,  we  develop  more  systematic  models  using  Proper  Orthogonal  De¬ 
composition  (POD)  and  Galerkin  projection.  These  models  are  very  accurate  for  the  con¬ 
ditions  they  are  calibrated  for,  but  they  do  not  generalize  well  to  other  angles  of  attack,  or 
more  complex  unsteady  maneuvers.  Thus,  while  the  models  in  this  section  are  interesting, 
we  expect  the  models  of  Section  [6]  to  be  more  useful  for  actual  implementation. 


4.1  Phenomenological  models 

Recall  from  Figure  [TO]  that  the  flow  over  a  wing  at  low  Reynolds  number  exhibits  a  Hopf 
bifurcation  at  a  critical  angle  of  attack,  at  which  the  flow  transitions  from  steady  separation 
to  unsteady  vortex  shedding.  The  simplest  system  of  differential  equations  that  describes 
this  behavior  is  given  by  the  normal  form  of  a  Hopf  bifurcation  |21|.  In  this  section,  we 
model  the  transient  and  steady-state  lift  associated  with  an  impulsively  started  2D  plate 
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using  this  simple  model,  along  with  a  decoupled  first-order  lag: 

x  =  (a  —  ac)nx  —  coy  —  ax(x2  +  y 2)  'j  r  =  r  [(ct  —  ac)y  —  ar2] 

y  =  (a-  ac)yy  +  lux  -  ay(x 2  +  y2)  >  =>•  0  =  lu  (22) 

i  =  —A  z  )  z  =  —Xz 

The  z  direction  is  decoupled  and  represents  the  exponential  decay  of  transient  lift  generated 

from  the  impulsive  start.  Transforming  the  (x,y)  system  into  polar  coordinates,  it  becomes 
clear  that  there  is  a  fixed  point  at  r  =  0.  This  fixed  point  undergoes  a  subcritical  Hopf 
bifurcation  at  a  =  ac  resulting  in  an  unstable  fixed  point  at  r  =  0  and  a  stable  limit  cycle 
with  radius  R  =  yj (a.  —  ac)fi/a.  The  limit  cycle  represents  the  fluctuations  in  lift  due  to 
periodic  vortex  shedding  of  a  plate  at  an  angle  of  attack  which  is  larger  than  the  stall  angle. 
Thus,  at  a  particular  angle  of  attack  cc,  the  unsteady  coefficient  of  lift  Cl  is  constructed 
from  the  average  lift  Cl  and  the  state  variables  y  and  z  as  follows 

Cl  =  Cl  +  y  +  z 


With  knowledge  of  the  actual  lift  vs.  time  from  numerical  experiment,  it  is  possible  to  tune 
the  parameters  (/u.,  A,  c u,a)  with  the  rates  of  decay,  period  of  shedding  and  amplitude  of 
stable  lift  fluctuations.  Initial  conditions  of  the  model  are  chosen  to  start  the  system  with 
the  right  transient  lift  and  phase.  By  properly  tuning  the  constants  with  experimental  data, 
this  model  will  closely  reproduce  the  transient  lift  dynamics  of  a  stationary  plate  for  a  wide 


range  of  angles  of  attack  a.  Figure  19  shows  a  typical  example,  for  a  =  35°. 


It  is  important  to  note  that  the  model  (|22j)  is  specifically  chosen  to  exhibit  a  supercritical 
Hopf  bifurcation  as  the  angle  of  attack  a  varies  through  ac.  Such  a  bifurcation  is  observed  in 
simulations  and  experiments  (Figure [l0)),  so  for  a  fixed  set  of  parameters,  the  model  (22 )  may 
be  tuned  to  match  the  observed  behavior  over  a  range  of  angle  of  attack  a.  Theoretically, 
one  would  expect  such  a  model  to  be  valid  only  for  values  of  a  close  to  the  bifurcation 
value  ac,  but  in  practice  we  have  observed  that,  even  for  fixed  parameters,  the  model  (22) 


may  be  tuned  to  match  simulations  reasonably  closely  over  a  wide  range  of  parameter  values, 
from  a  =  0  to  at  least  35°. 

It  also  is  possible  to  model  the  transient  lift  from  the  impulsive  start  by  coupling  z  and  (x,  y): 


x  =  (a  —  ac)yx  —  uiy  —  axz ' 
y  =  (a  —  ac)yy  +  lux  —  ayz 
ez=  -z  +  (x2  +  y2) 


r  =  r  [(a  —  afijy  —  az ] 

6  =  LU 


(23) 


ez  =  —  z  ■ 


For  e<  1,  trajectories  quickly  settle  to  the  “slow”  manifold  z  =  x2  +  y2,  which  reduces  to 
our  original  dynamics.  This  method  has  been  useful  in  characterizing  transients  in  the  wake 
of  a  cylinder |45|.  Because  the  flat  plate  at  high  angles  of  attack  is  a  bluff  body,  its  wake 
topology  should  be  structurally  equivalent  to  the  wake  behind  a  cylinder. 


The  models  (22)  and  (23)  capture  reasonable  dynamic  behavior,  but  they  have  a  severe 
limitation  for  the  present  application,  since  the  initial  values  of  the  states  (x,  y,  z)  need  to 
be  carefully  chosen  in  order  to  match  the  transient  lift  response.  In  practice,  this  transient 
response  is  precisely  what  we  are  interested  in  modeling,  but  these  initial  conditions  are 
not  known.  Instead,  these  states  are  excited  by  external  factors,  such  as  disturbances  or 
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Figure  19:  Coefficient  of  lift  vs.  time  for  stationary  plate  at  post-stall  angle  a  =  35°.  The 
solid  line  represents  the  oscillating  lift  curve  obtained  from  direct  numerical  simulation  and 
the  dashed  line  is  a  the  output  of  a  low  order  ODE  model  with  a  Hopf  bifurcation  and  a  lag. 
The  inset  panels  show  coherent  structures  in  the  flow  field  at  instances  of  high,  moderate 
and  low  lift. 
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changes  in  angle  of  attack  a  (which  can  lead  to  the  formation  of  leading-edge  vortices,  as 
stated  earlier).  One  idea  for  incorporating  these  effects  is  to  add  impulsive  forcing  to  the 
model  x  =  . . .  +  f(a,d),  so  that  the  terms  f  can  excite  the  states  x  through  a  change  of 
angle  of  attack.  However,  this  approach  was  not  pursued  in  this  project,  as  the  approaches 
of  Sections  [5]  and  [6]  seemed  more  fruitful. 

4.2  Models  using  POD  and  Galerkin  projection 

In  this  section,  we  summarize  our  results  determining  unsteady  aerodynamic  models  using 
Proper  Orthogonal  Decomposition  (POD)  and  Galerkin  projection.  These  should  not  be 
confused  with  the  models  obtained  using  Balanced  POD  (Section  [5j,  which  yielded  the 
superior  models  in  Section  [6] 

The  evolution  of  a  flow  field  via  direct  numerical  simulation  may  be  viewed  as  a  high  di¬ 
mensional  dynamical  system  u  =  X(u),  where  u  is  a  state  variable  representing  the  velocity 
components  at  each  spatial  location,  arranged  in  a  long  vector.  Therefore,  u  €  where 
N  is  the  number  of  grid  points  times  the  number  of  flow  variables. 

To  obtain  a  reduced  order  model,  it  is  first  necessary  to  construct  a  low  dimensional  sub- 
space  S  C  on  which  the  dynamics  may  be  projected.  Given  a  time  sequence  of  unsteady 
velocity  fields  {u^  from  DNS,  we  seek  a  projection  Pg  :  — >  S  so  that  the  pro¬ 

jection  error  jj  Y2h=i(\\uk  ~  Psuk ||)  is  minimized.  It  has  been  shown  that  this  minimization 
problem  is  equivalent  to  the  eigenvalue  problem 


Rip  =  Xip  where  R  =  XX*  and  X 


U\  U2  Um 


R  is  a  real,  symmetric  matrix  of  dimension  N  x  N  with  at  most  M  nonzero  eigenvalues. 
Therefore,  the  eigenvalue  problem  Rip  =  A  ip  is  equivalent  to  a  simpler  eigenvalue  problem 

Up  =  A  p  where  U  =  X*X 

of  dimension  M  x  M.  For  M  <C  N,  this  greatly  reduces  the  computation  and  is  known 
as  the  method  of  snapshots.  {pk}™=\  are  eigenfunctions  associated  with  the  rn  largest 
eigenvalues  of  R  (or  U)  and  are  known  as  POD  modes.  POD  modes  are  typically  computed 
after  subtracting  the  mean  flow  from  each  of  the  snapshots  u, .  Because  of  the  form  of 
the  minimization  problem  posed,  the  first  k  POD  modes  are  the  k  most  Energy-containing 
modes.  Although  energetic  modes  are  important,  it  has  been  shown  that  modes  including 
only  a  small  fraction  of  the  total  energy  can  be  dynamically  important  [26], 

Given  dynamics  u  =  X(u)  and  a  projection  Pg  onto  a  low  dimensional  subspace  S  C  K^, 
it  is  now  possible  to  project  the  discretized  Navier-Stokes  equations  onto  the  subspace  S, 
resulting  in  a  low  order  dynamical  system  model  for  the  full  equations  of  motion: 

m 

r  =  PsX(r),  where  r(t)  =  (p  +  ^  ak(t)pk 

k= 1 

The  dynamics  are  now  captured  as  a  low  dimensional  ODE  with  the  POD  mode  amplitudes 
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Mean  Flow 
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Figure  20:  Mean  flow  and  eigenvalues  of  POD  modes  for  stationary  flat  plate  at  a  =  45° 
and  Re=  100. 


as  variables: 


(r-X(r),tpk)  =  0 


(aj(pj,ifk)  -  (X(r),tpk)  =  0 


ak  =  (X(r),tpk)  for  k  =  l,...,m 


In  the  current  study,  ii  =  X(u)  is  the  momentum  equation: 

u  =  —{u  ■  X7)u  +  vV2u  —  Vp 
' - v - ' 

X(u) 


4-2.1  POD/Gcilerkin  model,  a  =  45° 

From  simulation  data  of  unsteady  flow  around  a  fixed  plate  at  Re  =  100  and  a  =  45°  which 
are  allowed  to  reach  steady  state  vortex  shedding,  POD  modes  are  computed  and  shown  in 
Figure  [2l]  In  this  configuration  the  plate  sheds  vorticity  periodically  from  the  leading  and 
trailing  edges.  By  first  subtracting  the  mean  flow,  it  is  possible  to  obtain  POD  modes  which 
are  in  energetic  pairs  as  seen  in  the  eigenvalue  plot,  Figure  [20]  Because  of  an  approximate 
convective  symmetry  in  the  periodic  shedding  case,  these  POD  modes  come  in  pairs  which 
appear  to  be  phase  shifted  by  7t/2. 


4-2.2  POD/Galerkin  model,  a  =  30° 


Figure  [23]  shows  POD  modes  for  the  unsteady  flow  around  a  fixed  plate  at  Re  =  100  and 
a  =  30°.  The  eigenvalue  plot,  Figure  22 
pairs  are  not  as  closely  matched. 


is  similar  to  the  a  =  45°  case,  except  that  the 


A  low-order  model  is  obtained  by  Galerkin  projecting  the  Navier-Stokes  equations  onto  these 
modes.  Using  the  projected  dynamics,  the  POD  mode  coefficients  are  integrated  forward  in 
time  and  used  to  reconstruct  an  approximate  flow  field.  From  the  approximate  velocity  field 
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Mode  1 


Mode  2 


Figure  21:  POD  modes  for  stationary  flat  plate  at  a  =  45°  and  Re=  100. 


Mean  Flow 


Eigenvalues 


Figure  22:  Mean  flow  and  eigenvalues  of  POD  modes  for  stationary  flat  plate  at  a  =  30° 
and  Re=  100. 
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Mode  1 


Mode  2 


Figure  23:  POD  modes  for  stationary  flat  plate  at  a  =  30°  and  Re=  100. 
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reconstruction,  an  FTLE  field  is  computed  which  agrees  very  well  with  the  FTLE  field  from 
DNS,  as  can  be  seen  in  Figure  [24]  The  comparison  of  POD  mode  amplitudes  between  DNS 


and  projected  dynamics  is  shown  in  Figure  25  It  is  interesting  that  the  model  agrees  well 
with  data  for  the  first  two  modes,  but  not  for  the  higher  order  modes.  However,  the  fact 
that  the  Lagrangian  coherent  structures  are  preserved  suggests  that  the  first  two  modes  are 
sufficient  for  reconstructing  an  accurate  model.  This  is  not  entirely  surprising  considering 
the  very  simple  sinusoidal  vortex  shedding  pattern  and  force  in  time.  For  higher  Reynolds 
number  flows,  the  lift  distribution  and  wake  structures  involve  higher  frequency  oscillations 
which  would  presumably  require  more  modes  to  approximate. 


Full  DNS 


Reconstructed 


Figure  24:  FTLE  field  for  DNS  vs.  reconstruction  from  Galerkin  projected  dynamics  for 
stationary  flat  plate  at  a  =  30°  and  Re=  100. 


5  Balanced  model  reduction 


The  models  developed  in  the  previous  section  demonstrate  good  agreement  for  the  conditions 
at  which  they  were  calibrated:  namely,  they  predict  the  unsteady  lift  forces  at  a  single  angle 
of  attack.  However,  the  phenomenological  models  of  section  4.1  are  obtained  in  an  ad  hoc 
fashion,  and  while  the  POD  modes  of  section  4.2  are  more  systematic,  it  is  not  clear  how  to 
extend  them  to  more  general  situations,  for  instance  in  which  one  has  a  gust  that  changes 
the  angle  of  attack  or  speed  of  the  freestream. 


In  order  to  represent  these  more  complex,  agile  maneuvers,  it  is  desirable  to  use  a  model 
reduction  procedure  that  includes  inputs  and  outputs :  for  instance,  disturbances  such  as  the 
instantaneous  angle  of  attack  or  freestream  speed  can  be  regarded  as  inputs,  and  desired 
quantities  such  as  lift  and  drag  forces  may  be  viewed  as  outputs.  Furthermore,  when  cast 
in  an  input-output  setting,  these  models  may  be  readily  used  for  control  design. 


Balanced  truncation  is  an  effective  technique  for  obtaining  reduced-order  models  of  input- 
output  systems.  However,  it  is  intractable  for  large  systems  such  as  fluids  problems.  An 
approximation  of  balanced  truncation,  called  Balanced  POD,  was  developed  by  Rowley  (49] . 
and  has  been  shown  to  produce  much  more  accurate  models  than  standard  POD/Galerkin 
models  [26] .  However,  a  criticism  of  this  technique  is  that  in  order  to  obtain  the  models,  one 
must  simulate  the  adjoint  equations ,  and  adjoint  simulations  are  not  always  readily  available. 
(Of  course,  they  are  never  available  for  experiments.) 
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Figure  25:  Mode  amplitudes  of  DNS  vs.  reconstruction  from  Galerkin  projected  dynamics 
for  stationary  flat  plate  at  a.  =  30°  and  Re=  100. 


In  the  following  subsection,  we  show  how  the  Balanced  POD  method  may  be  applied  even 
when  adjoint  information  is  not  available.  A  major  result  is  that  the  models  obtained  by 
Balanced  POD  are  precisely  those  obtained  by  an  existing  system  identification  method 
called  the  Eigensystem  Realization  Algorithm  (ERA). 

5.1  Balanced  POD  models  without  adjoints 

In  this  section,  we  summarize  the  steps  involved  in  approximate  balanced  truncation  (bal¬ 
anced  POD),  and  the  Eigensystem  Realization  Algorithm,  and  show  that  they  are  equivalent. 

Balanced  truncation  involves  first  constructing  a  a  coordinate  transformation  that  “balances” 
a  linear  input-output  system,  in  the  sense  that  certain  measures  of  controllability  and  ob¬ 
servability  (the  Gramian  matrices)  become  diagonal  and  identical  [43j.  A  reduced-order 
model  is  then  obtained  by  truncating  the  least  controllable  and  observable  states,  which 
correspond  to  the  smallest  diagonal  entries  in  the  transformed  system.  Unfortunately,  the 
exact  balanced  truncation  algorithm  is  not  tractable  for  the  large  state  dimensions  encoun¬ 
tered  in  fluid  mechanics.  However,  an  approximate,  snapshot-based  balanced  truncation 
algorithm,  referred  to  as  Balanced  Proper  Orthogonal  Decomposition  (balanced  POD)  was 
proposed  in  m,  and  has  been  used  successfully  in  several  examples  mmu- 
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The  second  technique,  the  eigensystem  realization  algorithm  (ERA),  has  been  used  both  for 
system  identification  and  for  model  reduction,  and  it  is  well  known  that  the  models  produced 
by  ERA  are  approximately  balanced  da  Esi.  Here  we  show  further  that,  theoretically, 
ERA  produces  exactly  the  same  reduced  order  models  as  balanced  POD.  This  equivalence 
indicates  that  ERA  can  be  regarded  as  an  approximate  balanced  truncation  method,  in  the 
sense  that,  before  truncation,  it  implicitly  realizes  a  coordinate  transformation  under  which 
a  pair  of  approximate  controllability  and  observability  Gramians  are  exactly  balanced.  This 
feature  distinguishes  ERA  from  other  model  reduction  methods  that  first  realize  truncations 
and  then  balance  the  reduced  order  models.  Note  that  in  ERA  the  Gramians,  and  the 
balancing  transformation  itself,  are  never  explicitly  calculated,  as  we  will  also  show  in  the 
following  discussions. 

For  both  techniques,  we  will  consider  a  high-dimensional,  stable,  discrete-time  linear  system, 
described  by 

x(k  +  1)  =  Ax(k)  +  Bu(k) 

y{k)  =  Cx(k), 

where  k  E  Z  is  the  time  step  index,  u(k)  E  Rp  denotes  a  vector  of  inputs  (for  instance, 
actuators  or  disturbances),  y[k )  E  M9  a  vector  of  outputs  (for  instance,  sensor  measurements, 
or  simply  quantities  that  one  wishes  to  model),  and  x(k)  E  denotes  the  state  variable 
(for  instance,  flow  variables  at  all  gridpoints  of  a  simulation).  These  equations  may  arise, 
for  instance,  by  discretizing  the  Navier-Stokes  equations  in  time  and  space,  and  linearizing 
about  a  steady  solution.  The  goal  is  to  obtain  an  approximate  model  that  captures  the  same 
relationship  between  inputs  u  and  outputs  y,  but  with  a  much  smaller  state  dimension: 


xr(k  +  1)  =  Arxr{k )  +  Bru(k ) 
y(k)  =  Crxr(k ) 

where  the  reduced  state  variable  xr(k)  E  Mr,  r  -C  n.  We  consider  the  discrete-time  setting, 
because  we  are  primarily  interested  in  discrete-time  data  from  simulations  or  experiments. 


5.1.1  Balanced  POD 


Here,  we  give  only  a  brief  overview  of  the  balanced  POD  algorithm,  and  for  details  of  the 
method,  we  refer  the  reader  to  [59].  The  algorithm  involves  three  main  steps: 


Step  1  :  Collect  snapshots.  Run  impulse-response  simulations  of  the  primal  system 
(24)  and  collect  mc  +  1  snapshots  of  states  x(k)  in  mcP  +  1  steps: 

X  =  [B  APB  A2PB  ■  ■  ■  AmcP B]  ,  (26) 


where  P  is  the  sampling  period.  In  addition,  run  impulse-response  simulations  for  the 
adjoint  system 

z(k  +  l)  =  A*z(k)  +  C*v{k)  (27) 

where  the  asterisk  *  stands  for  adjoint  of  a  matrix,  and  collect  m0  +  1  snapshots  of 
states  z{k)  in  m0P  +  1  steps: 


(A*)p  C* 


(A 


*\2 P  q* 


(A* 


\m0P 


c* 


(28) 


Calculate  the  generalized  Hankel  matrix, 


H  =  Y*X. 


(29) 
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Step  2:  Compute  modes.  Compute  the  singular  value  decomposition  of  H\ 


H  =  UYV* 


[Ui 


u2] 


Ei  O' 

0  0 

V* 

lV2  J 

UilhV? 


(30) 


where  the  diagonal  matrix  Si  E  MniXni  is  invertible  and  includes  all  non-zero  singular 
values  of  H ,  ni  =  rank(iL),  and  UfUi  =  V*Vi  =  Ini  XR] .  Choose  r  <  n±.  Let  Ur  and 
Vr  denote  the  sub-matrices  of  U \  and  V\  that  include  their  first  r  columns,  and  Er  the 
first  r  x  r  diagonal  block  of  Si.  Calculate 

d>r  =  XVrY.r'2]  d/r  =  YUrY~K  (31) 


where  the  columns  of  <Lr  and  dv  are  respectively  the  first  r  primal  and  adjoint  modes 

<r  • 


of  system  (24).  The  two  sets  of  modes  are  bi-orthogonal:  d/*$r  =  /rxr 


Step  3:  Project  dynamics.  The  system  matrices  in  the  reduced  order  model  (25)  are 
Ar  =  ^*A<5>r-  Br  =  d> *B ;  Cr  =  Cd»r.  (32) 


Note  that  the  n  x  n  controllability /observability  Gramians  are  approximated  by  the  matri¬ 
ces  XX*  and  YY* .  The  reduced-order  model  (25)  is  obtained  by  considering  a  subspace 
x  =  d>rxr,  and  projecting  the  dynamics  (24)  onto  this  subspace  using  the  adjoint  modes 
given  by  dV.  It  was  shown  in  [49]  that  d>r  and  'Lr  respectively  form  the  first  r  columns  of 
the  balancing  transformation/inverse  transformation  that  exactly  balance  the  approximate 
controllability/observability  Gramians  XX*  and  YY* . 


5.1.2  The  eigensystem  realization  algorithm 

The  eigensystem  realization  algorithm  (ERA)  was  proposed  in  |28]  as  a  system  identification 
and  model  reduction  technique  for  linear  systems.  The  algorithm  follows  three  main  steps 
[2811291: 


Step  1:  Run  impulse-response  simulations/experiments  of  the  system  (24)  for  (mc  + 
m0)P  +  2  steps,  where  mc  and  m0  respectively  reflect  how  much  effect  is  taken  for  con¬ 
sidering  controllability  and  observability,  and  P  again  is  the  sampling  period.  Collect 
the  snapshots  of  the  outputs  y  in  the  following  pattern: 


(CB,  CAB,  CAPB,  CAP+1B,  . . . 

\mcP  d  ri\rncP+lj^  (j j^(mc+m0)P ^  q ^{mc+rno)P+l  jjj'j 


CAmcPB,  CAr‘ 


(33) 


The  terms  CAkB  are  commonly  called  Markov  parameters.  Construct  a  generalized 
Hankel  matrix  H  E  M9(™o+i)xp(mc+i) 


H  = 


CB 

CAPB 


CAPB 

CA2PB 


CAm°pB  CA(m°+VpB 


CAmcPB 

CA(mc+l)P5 

(JA(rnc+m0)P 


(34) 


Step  2:  Compute  SVD  of  H,  exactly  as  in  (30),  to  obtain  Ui,  \j,  Sj.  Let  r  <  rank(iL). 
Let  Ur  and  Vr  denote  the  sub-matrices  of  U\  and  Vj  that  include  their  first  r  columns, 
and  Er  the  first  r  X  r  diagonal  block  of  Sj. 
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•  Step  3:  The  reduced  Ar,  Br  and  Cr  in  (25)  are  then  defined  as 


Ar  =  Er  ^UtH'VrY,r  2; 


1 

2  T/*  . 


Br  =  the  first  p  columns  of  S 2  V* ; 
Cr  =  the  first  q  rows  of  UrY 


(35) 


where 


H'  = 


CAB 


CAP+1B 


CAm°p+1B  CA<'m°+1)p+1B 


CAmcp+iB 


(JA{mc+m0)P+ 1 


(36) 


which  can  again  be  constructed  directly  from  the  collected  snapshots  (33). 


5.1.3  Theoretical  equivalence  between  ERA  and  balanced  POD 


The  first  observation  is  that,  with  X  and  Y  given  by  (26)  and  (]28|) ,  the  generalized  Hankel 


matrices  obtained  in  balanced  POD  and  ERA,  respectively  by  (]29[)  and  (34),  are  theoretically 


identical.  The  theoretical  equivalence  between  the  two  algorithms  then  follows  immediately: 


First,  H'  given  in  (36)  satisfies  H'  =  Y*AX,  which  implies  the  matrices  Ar  obtained  in  the 


two  algorithms  are  identical.  To  show  the  equivalence  of  Br,  first  note  that  the  SVD  (30) 

leads  to  E1  5 U(H  =  Y^V*,  which,  by  definition  of  Ur.  Vr,  Er,  implies  Yr  2 U*H  =  Yr V*. 
(Note  that  it  does  not  imply  H  =  Ur YrV* ,  since  UrU*  is  not  the  identity.)  Thus,  in  balanced 

POD,  Br  =  ^ZB  =  E 


r  2U*Y*B,  which  equals  the  first  p  columns  of  Er  2U*H  =  Y2  V* 


HVr Yr  2  =  UrY2 .  Thus,  in  balanced  POD, 

i  i 


i 

which  is  the  Br  given  by  ERA.  Similarly,  the  SVD  (30)  leads  to  HV iE1  2  =  U\Y(  and  then 


Cr  =  C&r  =  CXVrYir  2 ,  which  equals  the  first 


q  rows  of  HVr Er  2  =  UrY2 ,  the  Cr  given  by  ERA.  In  summary,  we  have: 

Main  result.  The  reduced  system  matrices  Ar.  Br  and  Cr  generated  in  balanced  POD  and 


ERA,  respectively  by  (32)  and  (35),  are  theoretically  identical. 

In  practice,  these  two  algorithms  may  generate  slightly  different  reduced  order  models,  be¬ 
cause  the  Hankel  matrices  calculated  in  the  two  algorithms  are  usually  not  exactly  the  same, 
due  to  small  numerical  inaccuracies  in  adjoint  simulations,  and/or  in  matrix  multiplications 
needed  to  compute  the  sub-blocks  in  the  Hankel  matrices.  In  the  following  discussions,  we 
compare  these  two  algorithms  in  more  detail. 

5.1-4  Comparison  between  ERA  and  balanced  POD 

While  ERA  and  balanced  POD  produce  theoretically  identical  reduced-order  models,  the 
techniques  differ  in  several  important  ways,  both  conceptually  and  computationally.  Neither 
ERA  nor  balanced  POD  calculate  Gramians  explicitly,  but  balanced  POD  does  construct 
approximate  controllability  and  observability  matrices  X  and  Y* ,  from  which  one  calculates 
the  generalized  Hankel  matrix  H  and  balancing  transformation.  Balanced  POD  thus  incurs 


additional  computational  cost,  because  one  needs  to  construct  the  adjoint  system  (27),  run 


adjoint  simulations  for  Y ,  and  then  calculate  each  block  of  H  by  matrix  multiplication. 
Thus  we  see  that  the  advantages  of  ERA  include: 
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Steps  in  computing 
reduced-order  models 

Approximate  time  (CPU  hours) 
balanced  POD  ERA 

1.  Linearized  impulse  response 

2 

4 

2.  Computation  of  POD  modes 

2 

2 

3.  Adjoint  impulse  responses 

30 

- 

(10  in  number) 

4.  Computation  of  the  Hankel  matrix 

7 

0.2 

5.  Singular  value  decomposition 

0.05 

0.05 

6.  Computation  of  modes 

1 

- 

7.  Computation  of  models 

0.02 

0.02 

Table  2:  Comparison  of  the  computational  times  required  for  various  steps  of  the  algorithms 
using  balanced  POD  and  ERA.  The  times  are  given  for  a  10-mode  output  projected  system. 
The  Hankel  matrix  is  constructed  using  (a)  200  state-snapshots  from  each  linearized  and 
adjoint  simulations  for  balanced  POD,  and  (b)  400  Markov  parameters  (outputs)  for  ERA. 


1.  Adjoint-free:  ERA  is  a  feasible  balanced  truncation  method  for  experiments,  since 
it  needs  only  the  output  measurements  from  the  response  to  an  impulsive  input.  Note 
that  ERA  has  been  successfully  applied  in  several  flow  control  experiments  [TTI  [TO] , 
as  a  system-identification  technique  rather  than  a  balanced-truncation  method.  In 
practice,  input-output  sensor  responses  are  often  collected  by  applying  a  broadband 
signal  to  the  inputs,  and  the  ARMARKOV  method  13137]  can  then  be  used  to  identify 
the  Markov  parameters,  or  even  directly  the  generalized  Hankel  matrix,  from  the  input- 
output  data  history. 


2.  Computational  efficiency:  For  large  problems,  typically  the  most  computationally 
expensive  component  of  computing  balanced  POD  is  constructing  the  generalized  Han¬ 
kel  matrix  H  in  (29),  as  this  involves  computing  inner  products  of  all  of  the  (large) 
primal  and  adjoint  snapshots  with  each  other.  ERA  is  significantly  more  efficient  at 
constructing  the  matrix  H  in  (|34|),  since  only  the  first  row  and  last  column  of  block 
matrices,  i.e. ,  the  (mc  +  m0  +  1)  Markov  parameters,  need  be  obtained  by  matrix 
multiplication.  All  the  other  mc  x  m0  block  matrices  in  H  are  copies  of  other  blocks, 
and  need  not  be  recomputed.  For  balanced  POD,  the  matrix  H  is  obtained  by  com¬ 
puting  all  the  ( mc  +  1)  X  ( m0  +  1)  matrix  multiplications  (inner  products)  between 
corresponding  blocks  in  Y*  and  X  in  (29).  Thus,  for  example,  if  mc  =  mo  =  200, 
the  computing  time  needed  for  constructing  H  in  ERA  will  be  about  only  1%  of  that 
in  balanced  POD.  See  Table  [3  for  a  detailed  comparison  on  computational  efficiency 
between  balanced  POD  and  ERA  in  an  example  of  the  flow  past  an  inclined  flat  plate, 
used  in  [39]. 


At  the  same  time,  balanced  POD  also  provides  its  own  advantages: 


1.  Sets  of  bi-orthogonal  primal/adjoint  modes:  Balanced  POD  provides  sets  of  bi- 
orthogonal  primal/adjoint  modes,  the  columns  of  <Lr  and  \Lr.  In  comparison,  without 
the  adjoint  system,  ERA  cannot  provide  the  primal  and  adjoint  modes.  At  best,  the 
primal  modes  may  be  computed,  using  the  first  equation  in  (31),  if  the  matrix  X  (26) 
is  stored  (in  addition  to  the  Markov  parameters).  But  the  adjoint  modes  cannot  be 
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computed  without  solutions  of  the  adjoint  system.  In  this  sense,  balanced  POD  incor¬ 
porates  more  of  the  physics  of  the  system  (the  two  sets  of  bi-orthogonal  modes),  while 
ERA  is  purely  based  on  input-output  data  of  the  system.  The  primal/adjoint  modes 
together  can  be  useful  for  system  analysis  and  controller/observer  design  purposes  in 
several  ways:  for  instance,  in  flow  control  applications,  a  large- amplitude  region  from 
the  most  observable  mode  (the  leading  adjoint  mode)  can  be  a  good  location  for  actu¬ 
ator  placement.  Also,  although  balanced  POD  is  a  linear  method,  a  nonlinear  system 
can  be  projected  onto  these  sets  of  modes  to  obtain  a  nonlinear  low-dimensional  model. 
For  instance,  the  transformation  x  =  <Lrxr,  xr  =  \I !*x  can  be  employed  to  reduce  a  full¬ 
dimensional  nonlinear  model  x  =  fix)  to  a  low-dimensional  system  xr  =  $*/($rrr). 
Finally,  if  parameters  (such  as  Reynolds  number  or  Mach  number)  are  present  in  the 
original  equations,  balanced  POD  can  retain  these  parameters  in  the  reduced-order 
models.  When  the  values  of  parameters  change,  the  reduced  order  model  by  balanced 
POD  may  still  be  valid  and  perform  well;  see  123  for  an  application  to  linearized 
channel  flow. 

2.  Unstable  systems:  Balanced  POD  has  been  extended  to  neutrally  stable  m  and 
unstable  systems  jT].  In  those  cases,  one  first  calculates  the  right /left  eigenvectors 
corresponding  to  the  neutral/unstable  eigenvalues  of  the  state-transition  matrix  A, 
using  direct/adjoint  simulations.  Using  these  eigenvectors,  the  system  is  projected 
onto  a  stable  subspace  and  then  balanced  truncation  is  realized  for  the  stable  subsys¬ 
tem.  ERA  for  unstable  systems  is  still  an  open  problem,  if  adjoint  operators  are  not 
available.  However,  we  note  that,  once  the  stable  subsystem  is  obtained,  ERA  can 
still  be  applied  to  it  and  efficiently  realize  its  approximate  balanced  truncation. 

ERA  for  systems  with  high-dimensional  outputs.  The  method  of  output  projection 
proposed  in  [49]  makes  it  computationally  feasible  to  realize  approximate  balanced  trunca¬ 
tion  for  systems  with  high-dimensional  outputs — for  instance,  if  one  wishes  to  model  the 
entire  state  x,  say  the  flow  held  in  the  entire  computational  or  experimental  domain.  This 
method  involves  projecting  the  outputs  onto  a  small  number  of  POD  modes,  determined 
from  snapshots  of  y  from  the  impulse-response  dataset.  This  method  can  be  directly  incor¬ 
porated  into  ERA  as  follows:  First,  run  impulse  response  simulations  of  the  original  system 
and  collect  Markov  parameters  as  usual.  Then,  compute  the  leading  POD  modes  of  the 
dataset  of  Markov  parameters  and  stack  them  as  columns  of  a  matrix  0.  Left  multiply 
those  Markov  parameters  by  0*  to  project  the  outputs  onto  these  POD  modes.  A  gener¬ 
alized  Hankel  matrix  is  then  constructed  using  these  modified  Markov  parameters,  and  the 
usual  steps  of  ERA  follow. 

5.2  Balanced  truncation  for  periodic  systems 

In  order  to  model  vortex  shedding  phenomena,  as  arise  for  wings  at  high  angles  of  attack, 
the  techniques  for  balanced  model  reduction  need  to  be  extended  to  periodic  systems.  This 
extension  is  straightforward,  but  somewhat  complex,  and  is  the  subject  of  the  current  section. 
In  this  context,  it  is  most  convenient  to  consider  discrete-time  systems,  which  may  be  viewed 
as  a  temporal  discretization  of  the  Navier-Stokes  equations. 

In  particular,  we  consider  linear  discrete-time  periodic  systems  of  the  form 

x(k  +  1)  =  A(k)x(k)  +  B(k)u(k);  y(k)  =  C(k)x(k),  (37) 
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with  state  x  E  Cn,  input  u  E  Cp,  output  y  E  C3,  and  T-periodic  matrix  coefficients 
A(-),  B(-) ,  C(-).  The  transition  matrix  in  (37)  is  :=  A(j  —  1 ) ^4(j  —  2)  •  •  •  A(i )  for  j  >  i, 
where  =  Inxn ■  Periodicity  implies  that  the  eigenvalues  of  Ftj+T,j)  are  independent  of  j. 

The  neutrally  stable  case  where  the  spectral  radius  p{F^j+x,j))  =  1  will  be  discussed  later. 
For  now,  assume  the  system  is  exponentially  stable ,  i.e.  p{F(j+T.j))  <  1-  The  controllability 
and  observability  Gramians  of  (37)  are  then  well  defined  and  are  T-periodic  in  j  [50]: 


3- 1 


wc(j)  ■■=  F(3,i+i)B^BWFtxi+iy 

i=— oo 
oo 

Wo{j) 


(38) 


i=j 


where  *  denotes  the  adjoint  operator. 

A  standard  lifting  procedure  [42]  recasts  (37)  in  T  input-output  (I/O)  equivalent  LTI  forms: 


Xj(t  +  1)  =  AjXj(t)  +  Bjujft)-, 
Vj(t)  =  CjXj(t)  +  DjUj(t ), 


(39) 


with  j  =  1  where  t  is  the  time  variable,  j  parameterizes  the  lifted  systems,  the 


state  Xj(t)  =  x(j  +  tT)  is  periodically  sampled  from  (37),  the  original  inputs  and  outputs 


T—l 


over  each  period  are  arranged  as  CpT  and  CqT  column  vectors  Uj(t)  =  [u(j  +  tT  +  i)]i=l 


and  yj(t)  =  [y(j  +  tT  +  i)]^1,  and  the  definitions  of  the  constant  matrices  Aj.  Bj ,  Cj  and 
Dj  readily  follow  from  the  variations  of  parameters  formula  in  (37),  e.g.,  Aj  = 
Assuming  exponential  stability,  the  controllability  and  observability  Gramians  of  the  j-th 
lifted  LTI  system  are 


i= 0 


F,o-.=  J2  (ty'c&A). 

i=0 


(40) 


The  following  statement  follows  from  the  periodicity  of  (37). 


Proposition  5.1.  Wjc  =  Wc(j)  and  WjQ  =  W0(j)  for  all  j  =  1, . . . ,  T. 


Proposition  5.1  enables  us  to  enjoy  the  best  of  both  worlds:  Whereas  lifting  enables  an 


appeal  to  LTI  balanced  truncation  in  the  lifted  domain,  as  discussed  through  the  remainder 
of  the  paper,  Gramian  computations  can  be  carried  in  the  original  periodic  setting,  where 
the  dimensions  of  the  input  and  output  spaces  are  much  lower:  p  and  q  instead  of  Tp  and 
Tq. 

5.2.1  Factorization  of  empirical  Gramians  using  snapshot-based  matrices 


the  exact  Gramians  are  substituted  by  approximate 


In  snapshot-based  methods 
empirical  Gramians  where  the  infinite  series  in  (38)  are  truncated  mm  [55]  at  a  finite 
m  <  oo: 

j- 1 


Wce(j-m):=  F(j,i+DB(i)BW*FU,i+ 1)! 


i=j—m 

j+m—l 


(41) 


Woe{j-,m)  :=  >  F, 

i=j 
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When  the  system  is  exponentially  stable,  truncation  is  justified  by  an  induced  norm  bound 
on  the  truncation  error,  obtained  by  a  geometric  series  argument  and  an  appeal  to  Propo¬ 
sition  15.11 


Lemma  5.2.  Assume  that  the  linear  periodic  system  (31)  is  exponentially  stable  and  let  m 
be  an  integer  multiple  of  the  period,  m  =  IT.  Then  the  following  induced  norm  error  bounds 
hold: 


\\Wc{j)-Wce{j-m)\\  ,  2 

HWcOOH  {j+T’j)l] 

\\W0(j)-W0e(j-,m)\\  ,  2 

\\W0(j)\\  (J+TJ)I1 


(42) 


Empirical  Gramians  can  be  factorized  using  snapshot-based  matrices. 

Proposition  5.3.  Let  B^\  i  =  1, . . .  ,p,  denote  the  i-th  column  of  B,  and  let  jW  ^  Cnxm 
be  defined  as 


F{jJ_m+2)B^(j  -  m  +  1), . . . ,  -  1) 


for  each  j  =  1, . . . ,  T  and  the  horizon  length  m.  Finally,  define  the  matrix  of  snapshots 


X (j ;  m)  :=  \x (1\j;m), . . . ,  X(p\j-?n) 


^  ^ nxmp 


Then  Wce(j;m)  =  X(j;m)X(j;m)* . 


As  illustrated  in  Figure  [26)(a) ,  the  columns  of  X (j ;  rn)  are  snapshots  of  impulse-response 
simulations  of  the  system  (|37|,  justifying  the  term  empirical  Gramian:  Invoking  the  T- 
periodicity  of  B(-)  and  F(-,-)(e.g.  ^(jj-m+T+i)  =  pQ-Tj-m+i));  one  observes  that  the 
m  columns  of  X^\m;j)  are  samples  at  times  j  —  kT,  k  =  0, ...,/  —  1  of  trajectories  of 
simulations  initiated  at  x(j  —  m  +  t)  =  B W  (j  —  rn  +  t  —  1) ,  t  =  1 , . . . ,  T,  assuming  m  =  IT. 
In  total,  Tp  simulations  and  mp  snapshots  are  needed  to  construct  X (j ;  rn) . 

An  analogous  observation  applies  to  the  empirical  observability  Gramian. 


Proposition  5.4.  Let  C^l\  i  =  1, . . . ,  q,  denote  the  i-th  row  of  C,  and  let  pW  e  Cnxm 
defined  as 


Y®(j\  m)  :=  [Flj+m_hj)C^(j  +  m-  1)*, 

F{j+m-2J)CU(j  +  m-2)*,...,Ci(jy 
for  each  j  =  1, . . . ,  T  and  the  horizon  length  m ..  Finally,  let 


be 


Y  (j ;  m) 


y(1)(j;m), . .  •  ,P(9)0';m) 


^  ^ nXmq 


Then  Woe(j;  m)  =  Y(j;  m)Y(j;m)*. 
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ic  =  [B(j  -  m)]W 

• — •- 


[B(j -m  +  T  - 1)]^ 


• — ► 


time  =  t(j  —  m  +  1)  t(j  —  m  +  T) 


(a) 


ic=  [C(j  +  m  —  l)*]w 

4 — • — ► 


ic=  [C(j  +  m  -  T)*]  W _ 


time  =  t(j  +  1) 


t(j  +  T) 


(b) 


•  m 


►  impulse-response 
simulations 


t(j-l)  i(i) 


• 


T  adjoint 
>  impulse-response 
simulations 


•  M 


t(j  +  m  -  1)  t(j  +  m) 


Figure  26:  (a)  The  T  impulse-response  simulations  corresponding  to  the  i-th  control  input, 
(b)  The  T  adjoint  impulse-response  simulations  corresponding  to  the  i-th  adjoint  control 
input. 


As  illustrated  in  Figure  26  T),  Y (j,  rri)  can  be  obtained  from  simulations  of  the  adjoint 
system 


z(k  +  1)  =  A{k)z{k)  +  C(k)v(k)  (43) 

where  k  =  j, . . . ,  j  +  m  —  1,  z  £  Cn,  v  £  Cq,  A{k)  :=  A(2j  +  m  —  k  —  1)*  and  C{k)  := 
C(2j  +  m  —  k  —  1)*.  By  periodicity,  Tq  adjoint  simulations,  and  in  total  mq  snapshots 
taken  at  time  j  +  kT,  k  =  1, . . . ,  l  are  needed  to  construct  Y  (j ;  m). 


5.2.2  Balanced  truncation  using  the  method  of  snapshots 

Fix  a  time  1  ^  j  ^  T.  Justified  by  Propositions |5.3|and|5.4| and  Lemma  5.2  let  X (j ;  mc )  and 
Y (j ;  m0)  be  computed,  allowing  mc  ^  m0,  as  factors  of  the  empirical  Gramians  Wce(j;  mc ),  Woe{j ;  mQ). 
By  Proposition  |5.1|  they  can  be  also  used  as  factors  of  the  empirical  Gramians  of  the 
j-th  lifted  system  (39).  The  method  of  snapshots  presented  in  |49|  then  leads  to  ap¬ 
proximate  balanced  truncations  in  the  lifted  LTI  setting,  as  follows:  Compute  the  SVD 
Y (j ;  m0)* X (j ;  mc)  =  UT,V*,  and  the  transformations  <L,  \L  that  exactly  balance  the  empir¬ 
ical  Gramians  of  the  lifted  system 


T  =  X(j-  mc) US”1/2;  ^  =  Y(j:  m0)U S'1/2. 


(44) 


Let  <Lr,  41,.  be  the  first  r  columns  of  T  and  \L,  comprising  the  leading  bi-orthogonal  balancing 
and  adjoint  modes  of  the  j-th  lifted  system.  (  Note  that  to  simplify  notation,  the  dependence 
of  <L,  \L,  <Lr,  \Lr  on  j  is  suppressed.  )  The  reduced  state  zj(t)  £  Cr  is  defined  by  the 
projection  Zj(t)  =  'L *Xj(t )  =  +  tT)  and  the  estimated  full  state  x(j  +  tT )  ~  <hr£,-(t). 
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The  reduced  model  of  order  r,  in  the  lifted  setting,  reads 


Zj(t  +  1)  =  'S>*AjQrZj(t)  + 

Vj(t )  =  Cj®rZj(t)  +  DjUj(t). 


(45) 


I/O  equivalence  of  (37)  to  the  lifted  (39)  means  that  the  reduced-order  system  provides  the 
sought  I/O  approximation  of  (37).  Note  that  improved  numerical  stability  of  the  computa¬ 
tions  above  can  be  achieved  by  first  representing  each  of  the  factors  X (j ;  ?nc)  and  Y (j ;  m0) 
in  terms  of  leading  orthogonal  bases,  obtained,  e.g.,  by  SVD  or  by  Krylov  methods. 


We  comment  in  closing  on  the  possibility  to  “un-lift”  the  reduced-order  lifted  system.  As 
discussed  in  |60|.  the  exact  Gramians  solve  an  allied  periodic  Lyapunov  equation,  thus 
providing  an  exact  periodic  balancing  and  an  “un-lifted”  balanced  truncation  in  the  periodic 
setting.  Using  the  method  of  snapshots,  There  are  two  computational  shortcoming  to  that 
approach  in  the  current  problem.  First,  the  computational  burden  is  high  when  T  3> 
1.  Second,  the  truncated  empirical  Gramians  used  here  do  not  form  an  exact  solution  of 
the  periodic  Lyapunov  equation.  Un-lifting  is  nonetheless  a  simple  task  if  the  balancing 
requirement  is  limited  to  the  periodically  sampled  system  (i.e. ,  to  a  lifted  system  for  one, 
fixed  j).  The  following  inductive  procedure  is  one  possible  solution:  Fix  4>(j)  :=  <l>r  and 
'L (j  +  T  —  1)  :=  \Lr.  Let  P(j  +  i)  be  the  rank-r  orthogonal  projection  on 
and  let  <I> (j  +  i  +  1)  =  \L(j  +  i)  G  Cnxr,  i  =  0, . . . ,  T  —  2,  satisfy  P(j  +  i)  =  <I> (j  +  i  + 
l)\L(jTi)*.  Then  a  periodic  realization  of  the  reduced  order  system  is  defined  with  Ar{k )  := 
tf(jfe)*A(lfe)$(fc),  Br(k)  :=  ty(k)*B(k)  and  Cr{k)  :=  C{k)${k). 

5.2.3  Output  projection  method 


The  computations  delineated  above  require  an  untenable  number  of  adjoint  simulations 
when  very  high  dimensional  outputs  are  considered;  e.g.,  when  the  output  is  set  identical 
to  the  state,  such  that  one  can  use  state  response  data  in  design  of  an  optimal  controller 
(e.g.  linear-quadratic  regulator)  or  to  analyze  system  dynamics  in  detail.  In  the  LTI  case 
[49j  proposed  to  project  the  output  on  the  (few)  leading  POD  modes  of  the  dataset  formed 
by  the  impulse  response  simulations.  Thus  one  invokes  the  kinematic  significance  of  POD 
modes,  to  reduce  the  dimension  of  the  output  space,  but  avoids  the  weakness  of  standard 
POD  models  that  use  them  as  dynamic  states.  Here  we  extend  the  output  projection  method 
to  periodic  systems. 


The  I/O  map  of  the  j-th  lifted  LTI  system  (39)  is  uctci  milieu  uy _ 

impulse-response  matrices  {Gj(t)}.  The  outp id- projected  lifted  system 


Xj{t  +  1)  =  AjXj(t )  +  Bjtij{t ); 

Vj (t)p  =  Pj  ( CjXj{t )  +  DjUj(t ))  , 


(46) 


is  designed  to  best  approximate  the  exact  impulse  response  of  the  original  lifted  system. 
Ideally,  the  low-rank  orthogonal  projection  matrix  Pj  should  thus  satisfy 


Pj  =  argmin 
1  Le%} 


J2\\Gj(t)  -  PjGjim' 


U= 0 


(47) 


where  Vf 


norm 


3  is  the  space  of  orthogonal  projections  of  rank  fop  <C  Tq.  When  the  Frobenius 
is  used  in  (47),  it  becomes  a  standard  projection  problem.  Its  solution  is 
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Pj  =  0j0*,  where  the  columns  of  Qj  are  the  leading  rop  POD  modes  of  the  datasets 

As  described  above,  the  optimal  Pj  is  generically  a  full  matrix.  Thus,  Vj(t)p  =  PjVj(t )  is  no 
longer  the  lifted  representation  of  the  output  of  a  periodic  system,  and  the  projected  system 
cannot  be  “un-lifted".  Rather,  for  each  t,  the  value  of  yj(t)p  is  determined  by  the  original 
response  along  an  entire  period.  In  particular,  we  lose  the  ability  to  compute  the  Gramian 


in  the  original  periodic  setting.  To  avoid  this  problem  we  impose  on  (47)  the  additional 
condition  that  the  projection  has  a  block  diagonal  form 


Pj  =  diag 


Pj(l),--- ,Pj(T) 


(48) 


where  each  q  X  q  diagonal  block  is  a  rank -rop  orthogonal  projection  with  rop  =  ropT.  This 
enables  to  un-lift  the  projected  lifted  system  I46h  to  an  output- projected  time-periodic  system 


(49) 


x{k  +  1)  =  A{k)x{k)  +  B(k)u(k ); 
y(k)P  =  P(k)C(k)x(k), 


where  the  T-periodic,  rank-rop  orthogonal  projection  P  is  defined  by  P(j  +  tT  +  i)  = 
P(j  +  i)  :=  Pj(i+  1),  i  =  0, . . .  ,T—  1.  The  constrained  optimization  problem  ( 47 )-( 48 )  is 


solved  as  an  equivalent  set  of  unconstrained  problems  in  the  periodic  setting,  invoking  the 
correspondence  oftheT,  q  x  pT  dimensional  blocks  of  Gj(t),  G(j+tT+i,j),  i  =  0,...,T— 1, 


to  the  impulse  response  of  (37),  as  detailed  in  |6j: 


Proposition  5.5.  Using  the  Frobenius  norm,  the  solution  of  the  constrained  optimization 


problem  (jl)  and  (fS)  is  equivalent  to  the  combined  solution  of  the  problems 


Pj{i  +  1)  =  _  argrnin  Y  \\G(j  +  tT  +  i,  j) 
{-P,(i+l)e'Prop}  \t=  o 

—Pj(i  +  1  )G(j  +  tT  +  i,j ) 


for  i  =  0, . . . ,  T  —  1. 


Proof.  By  a  reduction  to  a  standard  projection  problem. 


□ 


The  computation  of  the  structurally  constrained  optimal  Pj  of  the  form  (48)  is  thus  reduced 
to  T  unconstrained  optimization  problems  for  each  P(k),  k  =  j,  •  •  •  ,j+T—  1,  in  the  periodic 
setting.  Following  standard  POD  rationale,  the  solutions  are  P(k)  =  Q(k)0(k)* ,  where  the 
rop  columns  of  Q(k)  are  the  leading  POD  modes  of  the  dataset  {G(tT  +  0,  and  the 

approximation  error  between  the  output-projected  system  and  the  original  system  is 


J2\\Gj(t)  -  PjGj(t)\\2F 


t= o 

j+T—l  oo 

=  ^2  ^2  \\G(j  +  tT  +  i,j )  —  Pj{i  +  1  )G(j  +  tT  +  i,j ) 

i=j  t= o 
j+T-l  q 

=  J2  x^m 

i=j  m=rop+ 1 
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where  for  each  i,  A(z)i, . . . ,  A (i)q  are  the  descending-ordered  eigenvalues  G(tT+i,j)G(tT+ 
The  POD  modes  can  be  computed  by  the  method  of  snapshots  [56],  applied  to  datasets 
comprising  the  columns  of  the  impulse-response  matrices  {G(tT  +  i,j)}f= 0.  Conveniently, 
provided  that  mc  ^  (s+  1)T,  periodicity  implies  that  data  required  to  compute  these  snap¬ 
shots  have  already  been  obtained  during  the  computation  of  X (j;  mc),  as  described  in 


5.2.1 


For  instance,  the  matrix  C(j)X(j ;  mc)  includes  the  columns  of  matrices  {G(j  +  tT. 
The  empirical  factor  Y (j:  rn0)  of  the  corresponding  observability  Gramian 


Wop(j)  =  J2F{itj)C(i)*e(i)Q(i)*C(i)F{i!j) 


l=J 


is  needed  in  order  to  realize  the  snapshot-based  approximate  balanced  truncation  for  the 
output-projected  system  (49).  This  is  accomplished  with  only  Trop  (rop  -C  q)  impulse- 


response  simulations  of  the  adjoint  time-periodic  system  corresponding  to  the  output-projected 
system  (49),  whose  control  input  is  r op - d i me n s i o n al . 


In  closing  we  note  that,  for  additional  simplicity  and  a  requirement  of  a  single  SVD  com¬ 
putation,  one  can  also  use  a  single,  time-invariant  output  projection.  Under  this  constraint, 
the  optimal  selection  is  P  =  00*,  where  the  columns  of  0  are  the  leading  POD  modes 
of  the  entire  impulse-response  {{G(tT  +  1  °f  (37).  This  stronger  constraint 


implies  further  reduction  in  matching,  when  compared  with  the  optimal  solution  in  the  lifted 
domain. 

5.2.4  Summary:  procedures  of  balanced  POD  for  periodic  systems 


Following  the  terminology  in  |49],  the  approximate  balanced  truncation  method  for  linear, 
time-periodic  systems  is  termed  a  lifted  balanced  POD.  Its  main  steps  include: 


Step  0:  Fix  a  time  j,  1  j  ^  T,  as  the  time  point  for  lifting. 

Step  1:  Run  Tp  impulse-response  simulations  to  obtain  mcp  snapshots  and  form  the 


n  x  mcp  dimensional  X (j ;  rnc)  as  described  in  §  5.2.1 


Step  2:  Compute  y  =  Cx  from  stored  states  in  simulations  carried  to  compute 
X(j;mc).  Solve  for  the  POD  problems  for  the  periodically  sampled  y(j  +  tT  +  i), 
to  obtain  the  output-projection  matrices  0(j  +  z) ,  i  =  0,  ■  •  •  ,  T  —  1. 

Step  3:  Run  Trop  impulse-response  simulations  of  the  adjoint  output-proejcted  system, 


to  form  the  n  x  m0r0p  dimensional  matrix  Y (j :  rn0)  as  described  in  §  5.2.1 


Step  4:  Compute  the  SVD  of  Y (j:  rn0)*X(j:  mc )  and  the  balancing  modes  for  the  lifted 


system  given  by  (44). 


Step  5:  Compute  the  reduced  lifted  system  (45). 


Variants  include  skipping  Step  2,  when  the  output  dimension  q  is  small,  and  using  a  single, 
time-invariant  output  projection,  as  discussed  in  §  |5.2.3  The  reduced  system  can  be 
liftted  to  a  periodic  system,  e.g.,  as  described  in  closing 


5.2.2  As  in 


un- 

,  an  obvious 

dual  version  of  the  algorithm  addresses  the  case  of  a  high-dimensional  input  space,  with 
only  few  outputs.  This  case  is  motivated  by  systems  susceptible  to  distributed  disturbances, 
simultaneously  effecting  the  entire  state  (e.g.,  B  =  I). 
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5.2.5  The  neutrally  stable  case 


Consider  a  linear  periodic  system  (37)  that  arises  from  linearization  of  a  system  around  an 
asymptotically  stable  periodic  orbit.  By  Floquet  theory  |25|.  in  this  case  Aj  =  Fu+tj)  is  only 
neutrally  stable,  due  to  one  unity  eigenvalue  that  corresponds  to  persisting  perturbations 
along  the  periodic  orbit  in  the  linearization.  Balanced  truncation  cannot  be  directly  applied 
to  a  neutrally  stable  system,  as  the  infinite  series  used  to  define  Gramians  may  diverge. 

P]  presented  an  extended  version  of  balanced  POD  for  unstable  LTI  systems  that  have  small 
unstable  dimensions.  Following  the  idea  presented  in  [66j,  it  decomposes  the  system  dynam¬ 
ics  into  stable  and  unstable  parts.  Then  it  applies  approximate  balanced  truncation  to  the 
stable  dynamics  while  keeping  the  unstable  dynamics  exactly.  This  method  is  conceptually 
applied  here  to  periodic  systems  through  the  lifted  setting,  with  all  computations  executed 
in  the  periodic  setting.  First,  for  a  given  lifting  time  j,  define  a  projection  onto  the  stable 


subspace  Es  I  Aj  1  by  Pj  =  Inxn  —  jprp-  where  Wj,Vj  €  Cn  are  the  left/right  eigenvectors  of 


Aj  corresponding  to  the  unity  eigenvalue.  Dynamics  of  the  neutrally  stable  lifted  system  (39 ) 
is  thus  restricted  to  the  stable  subspace  of  Ar. 


Xj(t+l)s  =  AjXj(t)s  +  P jBjUj(t); 
Vj{t)s  =  Cj¥jXj{t)s  +  DjUj(t), 


(50) 


where  Xj(t)s  =  P jXj{t).  Lifted  balanced  POD  can  be  realized  to  this  projected  system 
describing  stable  dynamics.  Let  and  \L*  be  the  matrices  including  the  leading  rs 


balancing  and  adjoint  modes  of  the  projected  system  (50).  Then,  a  reduced  model  of  order 
•  s  +  1,  for  the  neutrally  stable  lifted  system  (39)  can  be  obtained  in  the  form  of 


r,  r  =  rs 


“I  T  ^  AjU  j 

(45),  where  now  <Fr  =  Vj\  ;  \Lr.  =  ^rs  w*Vj  ■  The  reduced  system  keeps  the  one¬ 

dimensional  neutrally  stable  dynamics  exactly,  while  the  exponentially  stable  dynamics  is 
reduced  to  the  order  of  rs. 

Numerically,  the  neutrally  stable  eigenvectors  of  Aj  can  be  calculated  using  a  Krylov  method, 


or  even  the  power  method:  By  running  a  control-free  simulation  of  the  periodic  system  (37) 
with  an  arbitrary  initial  condition  x^  ^  Es  (Aj) ,  one  can  approximate  Vj  by  x(j  +  IT) .  with 


a  large  l.  Similarly,  a  long-time  control-free  simulation  of  the  adjoint  periodic  system  (|43j)  is 
needed  to  approximate  Wj.  Then,  when  computing  the  transformations  and  for  the 


projected  system  (50),  one  follows  exactly  the  same  procedures  given  in  §  5.2.4  The  only 


5.2.1  the 


difference  is  that  in  the  Tp  simulations  of  the  periodic  system  (37)  described  in 
states  should  be  projected  onto  Es  (Aj^j  by  P j  at  time  j  —  m  +  T .  The  simulations  then 
resume  with  these  states  as  new  initial  conditions.  Similarly,  in  the  adjoint  simulations,  the 
adjoint  states  should  be  left-multiplied  by  P*  at  time  j  +  T  before  the  simulations  resume. 

By  construction,  this  method  is  applicable  to  other  neutrally  stable/unstable  periodic  sys¬ 
tems,  with  small  neutrally  stable/unstable  dimensions.  For  unstable  systems,  in  impulse- 
response  simulations  one  can  repeatedly  project  the  states  once  each  period,  using  P j,  to 
numerically  confine  the  dynamics  to  the  stable  invariant  subspace. 

5. 2. 6  Numerical  example 

To  illustrate  the  balanced  POD  algorithm,  consider  an  exponentially  stable  example  (sim¬ 


ilar  to  that  in  [14]):  a  linear  periodic  system  (37)  with  period  T  =  5,  state  dimension 
n  =  30,  output  dimension  q  =  30,  control  input  dimension  p  =  1,  and  {-A(fc)}|=1  are 
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randomly  generated  diagonal  matrices  with  diagonal  entries  bounded  in  [0.16, 0.96],  guaran¬ 
teeing  asymptotic  stability.  The  matrices  B(k)  and  C(k)  are  also  randomly  generated,  with 
entries  bounded  in  [0, 1]. 


Here  we  pick  the  “lifting  time”  j  =  1.  Choose  mc  =  m0  =  3 T  =  15.  Figure  [27|(a)  shows 
the  error  plots  of  the  infinity  norm,  ||G  —  GVHoo/HGHoo  versus  r,  the  order  of  the  reduced 
lifted  system.  Here  Gr  is  the  impulse-response  matrix  of  the  reduced  lifted  system  of  order 
r.  We  see  that  the  snapshot-based  balanced  truncation  gives  a  good  approximation  of  exact 
balanced  truncation.  Further,  the  balanced  POD,  even  with  low  orders  of  output  projection 
rop,  generates  satisfying  results.  Recall  that,  for  the  lifted  system,  the  order  of  output 
projection  is  fop  =  ropT. 


Figure  [27^b)  shows  comparisons  between  balanced  POD  results  with  the  same  order  of 
output  projection,  one  set  based  on  T-periodic  projection  matrices  along  one  period,  and 
the  other  using  single  time-invariant  projection  matrix  (see  §  5.2.3).  For  the  cases  where 


rop  are  low,  these  two  approaches  give  almost  identical  results,  or  even  the  latter  one  gives 
better  results.  However,  when  the  order  of  output  projection  rop  increases,  the  results  based 
on  T-periodic  projection  matrices  are  better  than  those  by  a  single  projection  matrix,  as  we 
expect. 


Figure  27:  Error  ||G  —  Gf Hoo/HGjloo,  for  lifted  balanced  POD  approach:  (a)  For  exact 
balanced  truncation(  ) ,  balanced  truncation  by  the  method  of  snapshots  but  without  output 
project  ion  (□) ,  balanced  POD  with  rop  =  1  (0),  balanced  POD  with  rop  =  3  (o),  and  the 
lower  bound  for  any  model  reduction  scheme  (— ).  All  output  projections  are  T-periodic. 
(b)  Time  varying  T-periodic  output  projections  versus  time-invariant  output  projections: 
balanced  POD  with  rop  =  1  (0),  balanced  POD  with  rop  =  3  (o)  and  balanced  POD  with 
rop  =  5  (+).  Solid  lines  correspond  to  cases  using  T-periodic  projection  matrices,  and  dashed 
lines  using  one  single  projection  matrix. 

This  algorithm  has  also  been  applied  to  a  neutrally  stable,  time-periodic  system  obtained  by 
linearizing  the  Ginzburg-Landau  partial  differential  equation  about  its  exponentially  stable 
time-periodic  solution;  see  m- 
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quasi-steady  &  added  mass 


U 


CL 


fast  dynamics 

Figure  28:  Schematic  for  reduced  order  model  of  Wagner’s  indicial  response. 


6  Control-oriented  extensions  of  classical  unsteady  models 

The  classical  models  described  in  Section[2]are  of  limited  use  for  control  design:  Theodorsen’s 
models  (Sec.  2.1)  are  valid  only  for  sinusoidal  maneuvers,  and  Wagner’s  models  (Sec.  2.2) 


require  a  cumbersome  convolution  integral. 

The  main  result  of  this  section  is  a  low-dimensional,  state-space  representation  of  unsteady 
aerodynamic  forces,  patterned  after  Wagner’s  indicial  response  model.  Step-response  simu¬ 
lations  in  either  angle  of  attack  a  or  vertical  position  h  are  performed  using  an  immersed 
boundary  projection  method.  Based  on  the  lift  coefficient  history  from  these  simulations,  the 


model  reduction  methods  discussed  in  Section  5.1  (specifically,  the  Eigensystem  Realization 


Algorithm)  is  used  to  construct  a  state-space  model  in  the  following  form: 


X 

Ar 

0 

0  ' 

X 

Br 

u 

= 

0 

1 

At 

u 

+ 

0 

u 

k-\- 1 

_  0 

0 

1 

u 

k 

At 

(51) 


CL(kAt)  =  [Cr  CLu  CL<1\ 


+  Cl iifc 


where  u  =  a  for  pitching,  u  =  h  for  plunging,  and  u  =  [a  h\ T  for  pitching  and  plunging. 

For  the  system  to  be  proper,  we  choose  ii  as  the  inpulQ  Furthermore,  u  and  u  appear 
explicitly  in  the  state,  since  there  is  a  quasi-steady  dependence  on  these  variables.  The 
additional  fast  dynamics,  modeled  by  G(s),  are  obtained  using  the  eigensystem  realization 
algorithm  (ERA),  resulting  in  a  low-dimensional  state-space  representation,  (Ar.  Br,  Cr),  of 


order  r.  This  is  shown  in  Figure  28 


The  output  of  the  model  is  the  lift  coefficient,  Cl,  and  the  model  structure  is  conveniently 
expressed  in  terms  of  the  stability  derivatives,  Clu.  •  Ui,  and  the  additional  fast  dynamics, 


*There  are  added  mass  forces  proportional  to  ii,  except  in  the  degenerate  case  of  pitching  about  the 
middle-chord. 
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Cr  ■  x: 


Cl  —  Clu  ■  u  +  Cl^  ■  u  +  Cxa  ■  u  +  Cr  ■  x 


(52) 


The  reduced-order  model  in  Eq.  (51 )  is  shown  to  agree  well  with  Wagner’s  indicial  response 
over  a  range  of  frequencies  and  maneuvers.  The  goal  of  this  section  is  not  to  demonstrate  the 
relevance  or  accuracy  of  Wagner’s  model,  but  rather  to  demonstrate  that  it  may  be  cast  in  the 
convenient  state-space  form,  either  for  use  itself,  or  as  the  foundation  for  a  more  sophisticated 
model.  However,  once  it  is  observed  that  the  reduced  order  model  faithfully  reconstructs 
Wagner’s  indicial  response,  comparisons  with  DNS  for  large-amplitude  maneuvers  provide 
interesting  context. 


6.1  Reduced-order  approximations  of  Wagner’s  indicial  response 

This  section  details  the  process  by  which  we  obtain  the  reduced  order  model  (ROM)  in  (51 ) 
for  the  lift  force  on  an  airfoil  as  a  function  of  motion  input  variables,  angle  of  attack  a 
and  vertical  position  h.  In  particular,  given  step-response  experiments,  we  first  subtract 


off  quasi-steady  and  added-mass  forces  in  Section  6.1.1  and  then  model  the  remainder  of 
the  fast  dynamics  using  the  eigensystem  realization  algorithm  (ERA)  in  Section  6.1.2  A 
step-by-step  summary  of  the  process  is  presented  in  Section  6.1.4[ 


Before  proceeding  with  the  details  of  the  method,  it  is  important  to  highlight  two  main 
characteristics  of  aerodynamic  step-response  experiments  that  motivate  the  form  of  the 
reduced  order  model.  First,  step  responses  in  certain  variables,  such  as  angle  of  attack  and 
vertical  velocity,  result  in  nonzero  steady-state  lift.  Second,  the  aerodynamic  response  to 
smoothed  step  inputs,  discussed  in  Section  |6.1.3  are  dominated  by  forces  proportional  to 
the  second  derivatives  of  the  configuration  variables,  a  and  h,  for  the  duration  of  the  step 
maneuver. 


The  case  of  pitching  about  the  middle  chord  is  an  exception  to  characteristic  2  above,  and 
a  step  in  vertical  position  violates  characteristic  1  above,  and  Section  |6.1.5  addresses  this 
special  case.  Finally,  the  theory  is  presented  for  the  case  of  single-input  single-output  (SISO) 
systems  in  one  configuration  variable,  either  a  or  h\  however,  the  theory  generalizes  to 


multiple-input  multiple-output  (MIMO)  systems,  which  is  briefly  discussed  in  Section  6.1.6 
6.1.1  Aerodynamic  impulse  response  and  stability  derivatives 


The  following  discussion  provides  a  theoretical  foundation  relating  the  various  transfer  func¬ 
tions,  state  space  realizations  and  impulse-response  simulations.  In  particular,  a  step  re¬ 
sponse  in  a  given  input  variable  u  may  be  viewed  as  an  impulse-response  in  the  time  deriva¬ 
tive  of  that  variable  u.  Moreover,  in  the  instance  when  the  output  y  depends  on  ii,  it  is 
necessary  for  the  input  to  be  ii,  since  there  is  no  way  to  represent  a  derivative  in  state-space 
form. 


In  the  discussion  that  follows,  it  may  be  helpful  to  think  of  u  as  angle  of  attack  a  and  the 
output  y  as  the  lift  coefficient  Cl-  However,  the  theory  is  presented  in  a  general  framework 
for  completeness. 

Impulse  response  in  u  According  to  linear  system  theory,  the  impulse-response  y^(t)  is 
the  output  corresponding  to  a  delta  function  input,  u (t)  =  5(t) jl]  By  linear  superposition,  it 

Mhe  output  is  the  lift  coefficient,  y  =  Cl  and  the  input  is  either  the  angle-of-attack,  u  =  q,  or  the 
vertical  position,  u  =  h. 
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Figure  29:  Transfer  function  Gu(s)  between  u  and  y. 


is  possible  to  construct  the  output  response  to  an  arbitrary  input  u (t)  by  convolution  with 
the  impulse-response  y^(t): 


y(t)=  f  yi(t -  rMT)dr  =  {y5u*u){t) 

Jo 


(53) 


Taking  the  Laplace  transform  of  both  sides  of  Eq.  ( 53 )  yields  the  convenient  form 


Y(s)  =  Gu(s)U(s) 


(54) 


where  Gu(s )  =  C  [t/u(t)]  is  the  transfer  function  relating  input  signals  U(s)  to  output  signals 
T(s).  This  is  shown  schematically  in  Figure  29 


Step  response  in  u  (impulse  response  in  u)  For  aerodynamic  systems,  an  impulse 
response  in  u  (either  angle  of  attack  or  vertical  position)  might  be  nonphysical.  Historically, 
it  was  practical  to  consider  the  step  response  y„(i)  (indicial  response)  to  step  inputs  u  = 
H(t),  where  H(t)  is  the  Heaviside  step  function.  The  step  response  in  u  is  the  same  as 
the  impulse  response  in  u,  so  that  =  y The  superposition  integral  in  Eq.  (53)  may  be 
rewritterjf] in  terms  of  the  step  response  ys(t ): 


y(t)  =  yS{t)u{ 0)  +  /  ys{t  -  t)u(t)cLt 

Jo 


(55) 


For  aerodynamic  systems,  there  may  be  a  nonzero  steady-state  value  for  the  step-response 
in  u.  The  steady-state  value,  given  by  Gu(0)  =  y^(oo)  corresponds  to  the  quasi-steady  lift 
slope  Clu.  It  is  convenient  to  split  the  step-response  y^(t)  =  y^(t)  into  the  steady-state  and 
transient  components: 

yl(t)  =  yi{  oo)  +  yi\t)  (56) 


so  that  Eq.  ( 55 )  becomes 


y(t)  =  yi(oo)u(0)  +  yi\t)u(0)  +  y^(oo)  /  u (t)cIt  +  I  y°J{t  -  r)u(r)dr  (57) 

Jo  Jo 

(58) 


,.8'r 


=  +  f  ySJ (t  -  r)u(r)dT 

Jo 

Again,  taking  the  Laplace  transform  of  both  sides  of  Eq.  (58)  yields 


*This  is  seen  by  substituting  y^t)  =  /()j/u(£  —  s)ds  into  the  following  integral  fg  y&(t  —  r)u(r)dT  and 
integrating  by  parts. 
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Figure  30:  Transfer  function  Gu(s)  between  u  and  y. 


Y(s)  =  Gu'{s)u(0)  +  yi(oo)U(s)  +  Gu'(s)£  [u(f)] 

=  +  sGu\s)  U{s ) 

where  Gf,(s)  =  C  [y£(t)]  is  the  transfer  function  from  u  to  y,  as  shown  in  Figure 
Rearranging  Eq.  (|60|)  yields  the  following: 


(59) 

(60) 
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Y(s)  = 


Gu(0) 


+  Gu'(s) 


sU(s) 


(61) 


GuW 


where  Gu(0)  =  y^(oo).  This  transfer  function  corresponds  to  the  following  state-space 
representation: 
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(62) 


where  (A,  B,C,  D)  is  a  state-space  realization  of  G^(s). 

Impulse  response  in  ii  It  is  observed  that  there  are  added-mass  forces  proportional  to 
ii.  These  added-mass  terms  appear  in  G^{s)  as  nonproper  derivative  terms,  so  that  there 
does  not  exist  a  state-space  realization  (A,  B,C,  D)  for  G^'(s).  Therefore,  we  must  choose 
ii  as  the  input  variable  for  the  system  to  be  proper  and  representable  in  state-space  form. 

After  subtracting  off  the  steady-state  lift  corresponding  to  Gu(0)  =  Clu,  it  is  possible  to 
integrate  the  impulse  response  in  u  to  obtain  the  impulse  response  in  ii: 


(63) 


This  amounts  to  multiplying  the  input  by  s  and  the  output  by  1/s.  It  is  now  possible 
to  repeat  the  process  above,  subtracting  off  the  steady-state  lift  Gu^O)  =  y-( oo),  now 
corresponding  to  Cl a.  The  result  is  a  transfer  function  Gu(s)  from  ii  to  y  as  shown  in 


Figure  31 
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Figure  31:  Transfer  function  Gu(s)  between  ii  and  y. 
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Let  (A,  B,C,  D)  be  a  state-space  realization  of  the  transfer  function  Gu'(s).  It  is  then 
possible  to  rewrite  the  time-domain  response  in  the  more  convenient  linear  system  framework 
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(65) 


T 

with  initial  condition  [0  u(0)]  .  It  is  worth  noting  that  the  x  dynamics  are  decoupled 

from  the  last  two  states  [u  u]  .  The  x  dynamics  model  the  transient  dynamics  from  the 
step-response,  while  the  additional  terms  Gu(0)u  and  account  for  the  quasi-steady 

lift  associated  with  constant  u  and  u.  The  next  section  details  an  approach  to  obtain 
an  r-dimensional  realization  (Ar,  Br,Cr,  Dr)  for  Ga'(s)  using  the  eigensystem  realization 
algorithm. 

6.1.2  Modeling  fast  dynamics  iLsing  the  Eigensystem  Realization  Algorithm 

One  method  of  obtaining  a  state-space  realization  for  the  transfer  function  Gy^s)  is  the 
eigensystem  realization  algorithm  (ERA).  As  described  in  Section |5.1[  this  method  involves 
taking  snapshots  (kAt)  from  an  impulse-response  function  y £  (t)  =  1[Gy,(s)]  and  it 

returns  a  reduced  order,  discrete-time,  state-space  model  with  time-step  At. 

Recall  that  we  begin  with  a  full-order  model,  and  wish  to  construct  a  reduced-order  model 
of  the  form 


x(fc  +  1)  =  Ax.(k)  +  Bu(k)  1 
y (k)  =  Cx(k )  +  Du(k)  J 


Reduction  Xr(/c  T  1)  —  ArXr(fc)  T  Br\l(k)  ,  , 

y  (k)  =  Cr-x.r{k)  +  Dru(k) 


Note  that  it  is  not  necessary  to  know  the  full-order  model  explicitly:  we  need  not  know  the 
matrices  (A,  B,  C,  D ),  but  we  do  need  to  be  able  to  obtain  impulse  response  data  from  the 
system.  This  data  may  be  obtained  either  from  simulations  or  from  experiments. 


The  steps  involved  are  described  in  more  detail  in  Section  5.1  but  here  we  give  a  brief 
summary: 
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1.  The  first  term  in  the  impulse  response  is  the  D  matrix,  or  feed-through  term.  If  the 
impulse  is  in  ii,  then  the  D  matrix  corresponds  to  the  term  Cxu  which  accounts  for 
added  mass  proportional  to  ii. 

2.  Gather  the  next  (mc  +  m0)  +  2  outputs  from  the  impulse-response  simulation,  where 
mc  and  m0  are  integers  representing  the  number  of  snapshots  for  controllability  and 
observability.  The  outputs  y (k)  =  CAk~lB  are  called  Markov  parameters,  and  are 
synthesized  into  a  generalized  Hankel  matrices: 


H  = 


CB 

CAB 

CAmcB 

CAB 

CA2B 

CAmc+1B 

CAm°B 

CAm°+1B  ■  ■  ■ 

(jAmc+m°B 

CAB 

CA2B 

CAmc+1B 

CA2B 

ca3b 

CAmc+2B 

(jAm°+1 

B  CAm°+2B 

(J  j^mc+m0+l  b 

H'  = 


3.  Compute  the  singular  value  decomposition  of  H : 

H  =  UZV*  =  [U{  U2\ 


Si  O' 

'V* 

0  0 

V* 
iV2  J 

=  U^Vf 


(67) 


(68) 


(69) 


4.  Finally,  let  Er  be  the  first  r  x  r  block  of  Si,  Ur,  Vr  the  first  r  columns  of  U\ ,  V\ .  and 
the  reduced  order  Ar,  Br ,  Cr  are  given  as  follows: 


X-^UfH'Vr^-112 

(70) 

first  p  columns  of  T.]J2Vf 

(71) 

first  q  rows  of  UrY}J2 

(72) 

6.1.3  Choice  of  step  function 


For  a  number  of  reasons,  an  actual  step  response  is  non-physical.  First,  it  is  impossible  to 
command  in  experiments  or  simulations,  because  it  would  correspond  to  a  body  instanta¬ 
neously  dematerializing  and  then  rematerializing  it  in  another  location.  An  alternative  is  to 
use  a  smoothed  step  maneuver  and  approach  the  limit  as  the  maneuver  becomes  very  rapid. 
As  the  maneuver  becomes  increasingly  rapid,  the  added-mass  forces  begin  to  dominate;  in 
fact,  a  good  rule  of  thumb  is  to  choose  a  maneuver  rapid  enough  that  the  lift  response  for 
the  duration  of  the  maneuver  is  dominated  by  added-mass. 

In  our  simulations,  the  duration  of  the  maneuver  is  T  =  .01  convective  time  units,  and 
the  amplitude  is  either  A  =  1°  ~  .01745  rad  in  the  case  of  pitching  or  A  =  .01745  chord 
lengths  in  the  case  of  vertical  position.  This  is  sufficiently  rapid  for  the  added-mass  forces 
to  dominate  during  the  maneuver. 


The  linear  ramp  maneuver  was  first  introduced  as  a  canonical  pitching  maneuver  to  compare 
and  study  various  experiments,  simulations  and  models jl6j.  The  equations  for  u  and  u  are 


given  in  Eqs.  (73-75),  and  the  maneuver  is  shown  in  Figure  32 
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Figure  32:  Linear  ramp  maneuver  in  u. 


Figure  33:  Diagram  illustrating  the  discrete-time  system  arising  from  pitch-ramp  step. 


G(t) 

u(t) 

u(t) 


=  log 


=  A 


cosh  (at) 


cosh(a(f  —  T)) 

G(t ) 


=  A 


max  G(t ) 

tanh(at)  —  tanh(a(t  —  T )) 
max  G(t) 


(73) 

(74) 

(75) 


The  start  of  the  maneuver  is  t  =  0  and  the  duration  of  the  ramp-up  is  T.  The  parameter 
a  effects  how  gradual  or  abrupt  the  ramp  acceleration  is.  By  choosing  a  large,  such  as 
a  =  1000,  it  is  possible  to  obtain  a  maneuver  where  the  u  acceleration  effects  are  localized 
near  time  t  =  0  and  t  =  T,  and  the  velocity  u  is  constant  throughout  much  of  the  maneuver. 
This  results  in  an  approximately  piecewise  linear  ramp-up,  where  the  transition  at  t  =  0 
and  t  =  T  are  smoothed. 

The  linear  ramp  maneuver  has  a  number  of  benefits  that  make  it  a  natural  choice  for  our 
smoothed  step  maneuver.  First,  the  boxy  profile  of  the  velocity  u  resembles  the  shape  of  a 
discrete-time  impulse  in  u  with  time  step  T.  This  makes  it  possible  to  run  simulations  that 
fully  resolve  the  maneuver  in  time,  and  then  sample  starting  at  the  middle  of  the  maneuver 
with  sampling  time  T  and  obtain  a  good  approximation  of  the  discrete  time  system.  This 
is  shown  schematically  in  Figure  [33| 
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6.1.4  Method  summary 
The  overall  method  is  summarized  as  follows: 

1.  Run  an  impulse-response  simulation  in  u,  which  is  a  step  response  in  u.  This  response 
exhibits  dominant  derivative  (ii)  term  during  the  smoothed  step  maneuver.  This 
indicates  that  the  input  should  be  ii  in  order  to  have  a  proper  system,  which  will  be 
addressed  in  step  3. 

2.  Subtract  off  CjJU  term,  corresponding  to  steady  state  lift  after  step  in  u: 

C'L  =  CL-  CLa  (76) 

3.  Integrate  to  get  impulse  in  il  (less  Clu  contribution) 

CL  =  J  C'Ldr  (77) 

4.  Subtract  off  Cl ,  term,  corresponding  to  steady  state  lift  due  to  constant  u, 

C'L  =  CL-  Cl,  (78) 

5.  Sample  smoothed  step  input  and  output  functions  starting  at  the  middle  of  the 
smoothed  impulse  with  sampling  time  At. 

6.  Obtain  Markov  parameters 

H  =  markov(Y, U , M) 

7.  First  term,  Hi  is  the  D  matrix  corresponding  to  il  feed-through  term,  Cl,  • 

8.  The  rest  of  the  Markov  parameters  Hj  go  into  ERA  code 

[Ar , Br , Cr , HSVs]  =  ERA(H,m,n,r) 

9.  Synthesize  into  final  form  of  ROM: 
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(79) 


The  preceding  algorithm  is  sufficiently  general  to  produce  reduced  order  models  for  an  airfoil 
pitching  about  various  points  along  the  chord,  plunging  vertically,  or  pitching  and  plunging. 
However,  because  of  the  exceptions  mentioned  in  the  introduction  to  this  section,  the  case 
of  middle-chord  pitching  and  vertical  plunging  require  minor  special  attention. 
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•  Exception  1.  Mid-chord  pitching:  Because  there  are  no  added-mass  forces  pro¬ 
portional  to  a  in  the  degenerate  case  of  pitching  about  the  middle  chord,  we  must 
manually  set  Cl &  =  0.  Although  the  value  obtained  in  the  algorithm  is  quite  small, 
any  non-zero  value  will  result  in  incorrect  behavior  at  high  frequencies  due  to  the 
relative  degree. 

•  Exception  2.  Vertical  plunging:  Given  a  step-up  in  vertical  position,  there  is  no 
corresponding  steady-state  lift,  and  so  we  manually  set  CLh  =0.  This  guarantees  the 
correct  behavior  at  low  frequencies. 


6.1.5  Possible  variations  for  plunging  and  pitching  about  mid-chord 


As  mentioned  above,  middle-chord  pitching  and  vertical  plunging  are  exceptions  and  require 
special  attention.  The  reduced  order  model,  Eq.  (51)  is  sufficiently  general  to  handle  these 
exceptions.  In  the  case  of  pitching  about  the  middle  chord,  we  manually  set  Cl&  =  0  and 
in  the  case  of  vertical  plunging  we  manually  set  CLh  =  0. 


Alternatively,  it  is  possible  to  eliminate  steps  3-4  above  for  mid-chord  pitching,  resulting  in 
the  model: 
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For  plunging,  one  may  use  a  step  in  vertical  velocity,  modifying  steps  1-4  and  results  in  the 
model: 
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(81) 


It  is  interesting  to  note  that  in  this  form,  with  the  input  for  pitching  about  the  mid-chord 
now  a,  these  models  should  be  equivalent  by  symmetry.  In  other  words,  a  plunging  motion 
h(t)  translates  directly  into  an  angle-of- attack  motion  aes(t)  by  a  change  of  variables. 

6.1.6  Multiple  Input  Multiple  Output  (MIMO)  generalization 

An  important  result  of  this  paper  is  that  the  above  algorithm  generalizes  to  MIMO  rep¬ 
resentations.  Lift  coefficient  dynamics  in  angle  of  attack  and  vertical  position  are  easily 
representable  in  a  combined  model. 

First,  treat  each  step  response  separately,  first  subtracting  off  Clu,  then  integrating,  sub¬ 
tracting  off  CLy .  and  finally  passing  through  a  Markov  parameter  identification  algorithm  to 
identify  Cl.,  and  the  transient  impulse  response  in  il.  Until  now,  we  have  exactly  repeated 
steps  1-7  individually  for  each  step  response.  Now,  the  ERA  algorithm,  step  8,  requires  a 
deviation.  The  Markov  parameters  CAkB  are  no  longer  coefficients,  but  matrices  of  size 


53 


Figure  34:  Step-responses  (left,  middle)  and  Hankel  singular  values  (right)  for  1°  pitch-up 
about  the  quarter  chord.  DNS  (black)  is  compared  with  a  6-mode  ERA  model  (red). 


input  x  output.  In  other  words,  we  stack  the  two  transient  impulse  responses  to  obtain  a 
multidimensional  impulse  response.  We  pass  this  through  the  ERA  algorithm  to  obtain  a 
representation  (Ar,  Br,Cr)  for  the  remainder  of  the  dynamics.  The  full  expression  for  the 
model  is: 
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CL(kAt)=[Cr  CLa  CLh  CL&  CLh\ 


6.2  Results:  small-amplitude  maneuvers 
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In  this  section,  we  present  the  results  of  the  procedure  described  in  section  6.1 
to  small-amplitude  pitching  and  plunging  maneuvers.  As  with  previous  examples  in  this 
report,  we  study  a  flat  plate  at  Re  =  100. 

6. 2. 1  Pitching 

Here,  we  present  results  for  a  flat  plate  pitching  about  the  quarter  chord.  The  results  are 
similar  for  airfoils  pitched  about  the  leading  edge  or  the  middle  chord,  with  the  exception 
that  for  middle-chord  pitching,  there  are  no  added  mass  effects  on  the  lift. 


Figure  34  shows  a  comparison  of  a  small-amplitude  step  for  the  full  simulation  (DNS), 


compared  with  a  6-mode  model  obtained  by  the  Eigensystem  Realization  Algorithm  (adjoint- 
free  Balanced  POD,  from  Section  [5|.  In  the  leftmost  plot,  the  responses  due  to  Ciaiphai 
Clci ,  and  CLa  have  been  subtracted  off,  as  described  in  Section  6.1.4  so  all  that  remains  is 
the  term  G(s )  in  Figure  28  The  middle  plot  shows  the  full  step  response  C'x(t),  with  these 


terms  included.  Note  that  in  both  cases,  the  agreement  of  the  6-mode  model  is  excellent. 
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Figure  35:  Frequency  response  of  Wagner’s  model  and  ERA  models  for  pitching  at  quarter 
chord. 


The  rapid  decay  of  the  Hankel  singular  values  in  the  rightmost  plot  shows  that  the  response 
may  be  captured  effectively  by  a  low-dimensional  model. 

Figure  [35]  shows  a  comparison  of  the  Bode  plots  (frequency  response)  for  the  full  simu¬ 
lation  (DNS),  compared  with  the  full  Wagner  model,  and  the  reduced-order  ERA  model. 
Theodorsen’s  model  is  also  shown,  for  comparison.  Note  that  in  order  to  compute  the  Bode 
plot  for  Wagner’s  model,  one  must  compute  a  convolution  integral  at  each  frequency,  while 
for  the  reduced-order  model,  one  has  a  simple  6th-order  transfer  function,  suitable  for  con¬ 
trol  design.  Agreement  of  this  reduced-order  model  with  Wagner  is  nearly  perfect,  and 
agreement  with  DNS  is  excellent  as  well. 


6.2.2  Plunging 


Here,  we  consider  a  maneuver  as  in  the  previous  section,  but  for  a  flat  plate  undergoing  pure 
plunge. 


Figure  [36]  shows  the  step  responses  for  the  plunge-up  maneuver,  comparing  the  full  DNS 
with  a  7th-order  ERA  model.  As  in  Figure  34  the  middle  plot  is  the  actual  lift,  while  the 
leftmost  plot  shows  the  lift  with  the  terms  due  to  CLdiphai  Clq,>  and  Clu  subtracted  off,  as 
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Figure  36:  Step-responses  (left,  middle)  and  Hankel  singular  values  (right)  for  plunge-up 
maneuver.  DNS  (black)  is  compared  with  a  7-mode  ERA  model  (red). 


described  in  Section  6.1.4 


The  step  responses  are  visually  identical. 


Figure  |37|  shows  the  corresponding  Bode  plots  for  these  cases,  including  a  comparison  with 
Wagner  and  Theodorsen’s  models.  All  of  the  models  match  the  DNS  well,  except  for  a 
15-20°  discrepancy  in  phase  for  very  low  frequencies.  As  before,  the  ERA  model  agrees 
nearly  perfectly  with  Wagner’s  full  indicial  response,  but  is  much  more  convenient  to  use  for 
control  design. 


6.3  Results:  large-amplitude  maneuvers 


Finally,  we  apply  these  models  to  a  canonical  unsteady  maneuver  considered  in  [46j.  In  this 
maneuver,  the  airfoil  pitches  up  to  a  large  angle  of  attack  of  45°,  holds  for  one  convective 
time  unit,  and  then  pitches  back  down  to  zero  angle  of  attack.  Wagner’s  models,  and  the 
approximations  to  them  given  in  the  previous  section,  are  designed  for  small  angles  of  attack, 
and  are  not  meant  to  capture  such  large  excursions.  However,  we  compare  the  model’s 
predictions  here  to  better  understand  the  significance  of  nonlinearities  on  the  unsteady 
aerodynamics. 


Figure  [39]  shows  the  lift  when  the  airfoil  is  pitched  about  the  quarter  chord,  at  Reynolds 
number  100.  The  lift  from  the  full  DNS  is  compared  with  the  predictions  of  Wagner’s  model, 
and  our  6-state  low-order  model  from  Section  6.2.1  The  agreement  between  our  model  and 
Wagner  is  nearly  perfect,  as  expected.  The  agreement  with  DNS  is  good  up  to  about  2 
convective  time  units,  at  which  time  the  angle  of  attack  exceeds  20°.  The  large  peaks  in  lift 
are  due  to  the  added  mass  response,  and  the  models  are  very  effective  at  capturing  these. 
The  agreement  is  poor  between  about  2  and  4  convective  time  units,  but  then  matches  quite 
well  as  the  airfoil  pitches  back  down  to  zero  angle  of  attack. 


One  of  the  reasons  for  the  discrepancy  for  large  angles  of  attack  is  that  in  the  linearized 
models,  the  force  is  normal  to  the  airfoil’s  surface,  while  the  lift  is  defined  as  the  force  normal 
to  the  freestream  velocity.  Thus,  a  simple  cos(a)  correction  can  be  applied  to  correct  for  this 
difference.  When  this  correction  is  applied  (shown  in  the  right  plot  in  Figure  39  agreement 
with  DNS  is  good  for  the  entire  pitch-up  maneuver,  but  then  begins  to  disgree  during  the 
hold  and  pitch-down  portions. 


F igur e |40| shows  the  response  analogous  to  that  shown  in  Figure[39[  but  for  the  airfoil  pitched 
about  the  middle  chord.  Here,  there  are  no  added  mass  forces  from  the  unsteady  pitching,  so 
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Figure  37:  Frequency  response  of  Wagner’s  model  and  ERA  models  for  plunging. 
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Figure  38:  Canonical  maneuver  with  pitch-up,  hold,  and  pitch-down  phases,  showing  the 
instantaneous  vorticity  field  at  various  points  along  the  maneuver,  at  Re  =  100. 


Figure  39:  Comparison  of  lift  model  for  canonical  pitch-ramp-hold  maneuver  about  quarter 
chord :  Left:  linearized  model  prediction;  Right:  model  prediction  with  cos(a)  correction. 
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Figure  40:  Comparison  of  lift  model  for  canonical  pitch-ramp-hold  maneuver  about  middle 
chord:  Left:  linearized  model  prediction;  Right:  model  prediction  with  cos(a)  correction. 


the  response  looks  significantly  smoother.  As  before,  the  linearized  model  prediction  is  good 
when  the  angle  of  attack  is  small,  but  deviates  when  the  angle  becomes  large.  The  cos(a) 
correction  improves  the  response  during  the  pitch-up  portion,  but  degrades  the  response 
during  the  pitch-down  portion. 


6.4  Conclusions 


The  goal  of  this  section  is  to  present  a  procedure  for  constructing  unsteady  aerodynamic 
models  that  are  suitable  for  control  design.  The  models  are  based  on  Wagner’s  indicial 
response  |63j ,  and  agree  with  Wagner’s  models  to  an  arbitrarily  high  degree  of  accuracy. 
While  Wagner’s  models  represent  the  lift  as  a  convolution  integral  that  is  cumbersome  to 
compute,  and  not  suitable  for  control  synthesis,  our  procedure  produces  state-space  models 
that  may  be  used  directly  for  control  design.  Furthermore,  our  models  are  formulated 
in  a  way  that  directly  builds  upon  standard  approaches  incorporating  classical  “stability 
derivatives”  Clu,  Cl&,  Cl&  (as  in  mi  as  shown  in  Figure  [28] 


Our  models  predict  the  unsteady  response  very  well  for  small-amplitude  maneuvers,  as 
shown  by  the  frequency  responses  in  Figures  [35]  and  37  For  large  amplitude  maneuvers, 


such  as  the  canonical  maneuver  shown  in  Figure  [38]  agreement  is  poor  when  the  angle  of 
attack  becomes  larger  than  about  20°.  More  work  is  needed  to  extend  these  models  to  the 
nonlinear  regime. 
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